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Abstract 

This work deals with the simulation of Wishart processes and affine diffusions on positive semidefinite 
matrices. To do so, we focus on the splitting of the infinitesimal generator, in order to use composition 
techniques as Ninomiya and Victoir [20] or Alfonsi pQ. Doing so, we have found a remarkable splitting 
for Wishart processes that enables us to sample exactly Wishart distributions, without any restriction 
on the parameters. It is related but extends existing exact simulation methods based on Bartlett's 
decomposition. Moreover, we can construct high-order discretization schemes for Wishart processes and 
second-order schemes for general affine diffusions. These schemes are in practice faster than the exact 
simulation to sample entire paths. Numerical results on their convergence are given. 
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Introduction and first definitions 

This paper focuses on simulation methods for Wishart processes and more generally for afhnc diffusions 
on positive semidefinite matrices. Before explaining our motivations and our main results, we start with a 
short introduction to these processes. Even though we use rather standard notations for matrices, they are 
recalled in Appendix [5] and we invite the reader to give first a quick look at it. Wishart processes have 
been initially introduced by Bru [3J 0] . They are also named because their marginal laws follow Wishart 
distributions. Very recently, Cuchiero et al. [5] have introduced a general framework for affine processes on 
positive semidefinite matrices Sj~([R) that embeds Wishart processes and includes possible jumps. In this 
paper, we only consider continuous processes of this kind. Such processes solve the following SDE: 



Xf=x+f (5 + B{Xf)) ds + f 
Jo Jo 



(1) 



Here, and throughout the paper, (Wt,t > 0) denotes a d-by-d square matrix made of independent standard 
Brownian motions, 

x,a € Sj(R), a G M d (R) and B e C{S d {R)) (2) 
is a linear mapping on iS<j(R). Wishart processes correspond to the case where 

3a > 0,a = aa T a and 3b e M d (R),yx G S d (R),B(x) = bx + xb T . (3) 
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When = 1 , ([T]) is simply the SDE of the Cox- Ingersoll- Ross process that has been broadly studied, and we 
will implicitly assume that d > 2 throughout the paper. Weak and strong uniqueness of SDE ([T]) has been 
studied by Bru g], Cuchiero et al. [B] and Mayerhofer et al. [15] ■ We sum up here their results. 

Theorem 1 — If x € <Sj"(R), a - (d - l)a T a e iSj~(R) and B satisfies the following condition 

Vxx,x 2 eS+(R), Ti( Xl x 2 ) = =>■ Tr(B( Xl )x 2 ) > 0, (4) 

there is a unique weak solution to the SDE ([TJ. We denote by AFFd(x,a,B,a) the law of (Xf) t >o and 
AFFd(x,a, B, a;t) the marginal law of Xf . If we assume moreover that a — (d + l)a T a G <S^(R) and 
x fc Sj'*(R), there is a unique strong solution to the SDE JTJ). 

Under the parametrization of Wishart processes ([3]), condition Q is satisfied and weak uniqueness holds 
as soon as a > d — 1. In that case, we denote by WISd(x,a,b,a) the law of the Wishart process (Xf) t >o 
and WISd(x,a,b,a;t) the law of Xf . 

Throughout the paper, when we use the notation AFFd(x, a, B, a) or AFFa(x, a, B, a; t) (resp. WISd{x, a, b, a) 
or WISd(x, a, b, a; t)), we implicitly assume that a— (d— l)a T a € <5>j"(R) (resp. a > d— 1) and B satishes flU 
so that weak uniqueness holds. 

In her Ph.D. thesis [3J, Bru has introduced Wishart processes and used them in biology to study perturbed 
experimental data. Recently, a great attention has been paid to Wishart processes for applications in finance. 
Namely, Gourieroux and Sufana [13] and Da Fonseca et al. [7] have suggested to use these processes to model 
the instantaneous covariance matrix of d assets. It naturally extends stochastic volatility models for only one 
asset like the Heston model [15j . Obviously, processes on positive semidefinite matrices are really interesting 
to model the evolution of a dependence structure because they can describe a covariance matrix. However, 
when dealing with applications, it is in general crucial to be able to sample paths of such processes and make 
Monte-Carlo algorithms. 

To the best of our knowledge, there is few literature on simulation methods for Wishart and general 
affine processes ([1]). Wishart distributions have been intensively studied in statistics when a G N. In this 
case, exact simulation methods have proposed by Odell and Feiveson [21], Smith and Hocking [22] and 
Gleser [TT] to mention a few. Concerning discretization schemes, the usual Euler-Maruyama scheme is not 
well-defined because of the square- root. This is what already happens for the Cox-Ingersoll-Ross process 
which corresponds to the case d = 1. One has then to find specific schemes. Recently, Benabid et al. [2] 
and Gauthier and Possamai' [9] have proposed numerical approximations for Wishart processes that are 
well defined under some restrictions on the parameters. However, there is no result on the accuracy of 
their methods. Currently, Teichmann |25] is working on dedicated schemes for general affine processes by 
approximating their characteristic functions. Our study here is only dedicated to the diffusion ([I]). 

Initially, our goal was to find high order discretization schemes for Wishart processes by splitting operators 
and using scheme compositions. Indeed, this approach has already proved to be very efficient for other affine 
diffusions, see pQ. Doing so, we incidentally have found a remarkable splitting for some canonical Wishart 
processes which enables us to sample exactly Wishart processes. In particular, our result extends the 
Bartlett's decomposition that is commonly used to sample central Wishart distributions when a E IN. This 
splitting also enables us to get high order discretization schemes for Wishart processes. Then, by using 
scheme composition, we also get a second order scheme for any affine diffusion ([1]). Contrary to the exact 
scheme, one has then to study the discretization error, and we provide a rigorous analysis of the weak error. 

This paper is structured as follows. First, we present some general results on affine diffusions. We 
calculate their infinitesimal generator and obtain interesting identities in law that are intensively used next 
for the different simulation methods. Section [2] is devoted to the exact simulation of Wishart processes. It 
exhibits a remarkable splitting of the infinitesimal generator and shows how it can be used to sample exactly 
any Wishart distribution. Section [5J deals with high order schemes for affine diffusions. Thanks to the 
remarkable splitting, we are able to construct a third order scheme for Wishart processes and second order 
schemes for affine diffusions. Last, we give numerical illustrations of our convergence results in Section 0] 
We compare the time required by each method and also give a possible application of our results in finance. 
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1 General properties of affine processes on positive semidefinite 
matrices 



1.1 The infinitesimal generator on Aid(R) and <Sd(IR) 

We start this subsection with a simple Lemma, which is useful to calculate the infinitesimal generator of 
affine processes. Its proof is basic and is left in Appendix IB. II 

Lemma 2 — Let (J r t)t>o denote the filtration generated by (Wt,t > 0). We consider a process (Yt)t>o valued 
in G 5(j(R), and we assume that there exist continuous (J- t )-adapted processes {A t )t>o, {Bt)t>o and {Ct)t>o 
that are respectively valued in M.d{R), A^d(R) and <S<j(R) so that Y t admits the following semimartingale 
decomposition: 

dY t = C t dt + B t dW t A t + AjdW^Bj (5) 
Then, for i,j, m, n 6 {1, . . . , d}, the quadratic covariation of (Yt)ij and (Yt) m ,n is given by: 

d{{Y t )i : j, {Y t ) m . n ) = (B t Bf ) i _ rn (AjA t )j. n + (B t Bj') i _ n (ATA t ) Jtm + (B t Bf)j_ rn (AjA t )^ n +(B t Bj)j_ n 

(6) 

It is worth to notice that the quadratic covariation given by ([5]) depends on A t and B t only through the 
matrices AfA t and B t Bj . Lemma [5] enables us to calculate easily the infinitesimal generator for the affine 
process (QJ which is defined by: 

x e S+(R), L M f{x) = Urn SZiS ~ for f e C 2 {M d (R), R) with bounded derivatives. 



Proposition 3 — Infinitesimal Generator on Alrf(R). Let (X^) t >a ~ AFFd(x, a, B, a) be an affine 
process. Its infinitesimal generator on VWd(R) is given by: 

L M = Tr([a + B{x)]D M ) + ^{2Ti(xD M a T aD M ) + Tv{x(D M ) T a T aD M ) + Tr{xD M a T a(D M ) T )}, (7) 

where D M = (9i,j)i<i,i<d- 

Proof : On the one hand, the drift part of the operator is given by l=1 ctk,idk,i J r^'k i=i(^( x ))k,idk,i = 
Tr((a + B(x))D M ). On the other hand, we get from ([6]) that the diffusion part of the operator is: 

+ 2 E»,m,i.pl x j,m{ a a)i,ndi,jdm,n + 2 En.m,i. j ^j." ( a a )i,m.dijd m ^ n 

= \{Tx{xD M a T a{D M ) T ) + Ti{x{D M ) T a T a{D M ) f ) + Tr(xD M a T aD M ) + Tr (x (D M ) T a T aD M )} 
= l[Tr(x(D M ) T a T a{D M ) T ) + 2Tr(xD M a T aD M ) + Tr(x{D M ) T a T aD M )}. 



Here, we have given the infinitesimal generator on J\4 d (R), while we know that the affine process (Xf )t>o 
takes values in 5+(R) C S d (R). Thus, we can also look at the infinitesimal generator of this diffusion on 
Sd(R), which is defined by: 

x e Sj(E), L s f{x) = lim E ^^ Xt U - IM f or / e C 2 (5 d (R), R) with bounded derivatives. 
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For x G <Sd([R), we denote by xuj\ — Xij — Xj t i the value of the coordinates and so that 

x = Y^i<i<j<d x {%j}( e d :> + ^i^i^d ) ( see notations in Appendix |A"| . For / € C 2 (S d (R), R), we then denote 
by d{ij}f its derivative with respect to the coordinates xuj\. We introduce 

tt: M d (R)~>S d (R) 
x n- (x + x T )/2, 

that is such that tt(x) = x for x € Sd(R)- Obviously, / o 7r G C 2 (.Md(R), R) and we have 

L s f(x) = L M f 0ir (x). 

By the chain rule, we have for a; G <Sd(R), dijfoTr(x) = (l,=j + ^ijtj)9uj\f(x) and get the following result. 



Corollary 4 — Infinitesimal Generator on <S(j(R). TTie infinitesimal generator on S d (R) associated to 
AFF d (x,a, B,a) is given by: 

L s = Tr([a + J B(x)]L>' s ) + 2Tr(>i7 s a T aD s ), (8) 
where D s is defined by Df j = (li=j + , for 1 < i,j < d. 



Of course, the generators L M and L s arc equivalent: one can be deduced from the other. However, L s 
already embeds the fact that the process lies in <5><j(R), which reduces the dimension from d 2 to d(d + l)/2 
and gives in practice shorter formulas. This is why we will mostly work in the sequel with infinitesimal 
generators on ^(R). Unless it is necessary to make the distinction with L , we will simply denote L = L s . 



1.2 The characteristic function of Wishart processes 

As for other affine processes, the characteristic function of affine processes on positive semidefinite matri- 
ces can be obtained by solving two ODEs. In the case of Wishart processes, it is possible to solve explicitly 
these ODEs by solving a matrix Riccati equation (see Levin [E]). Here, we give the closed formula for the 
Laplace transform and a precise description of its set of convergence. 



Proposition 5 — Let Xf ~ WISd{x, a, b, a; t), q t = L exp(sb)a T aexp(sb T )ds and m t — exp(tb). We 
introduce the set of convergence of the Laplace transform of Xf , T>b ta . t = {u G S d (R), E[exp(Tr(wATf ))] < oo}. 
This is a convex open set that is given explicitly by 

% ; t = {^5d(R),Vse [0,t],I d -2q s veg d (R)}. (9) 

Besides, the Laplace transform of Xf is well-defined for v = vr + ivj with vr G T>b, a - t ,vi G <S<z(R) and is 
given by: 

E[exp(Tr(^f ))] = ^™f< - 2 qt v^rn t x m J]) 

aet(I d - 2q t v) 2 

The characteristic function corresponds to the case vr = that clearly belongs to 2\ a; t- The proof of this 
result is given in Appendix IB. 21 Let us remark here that for Xf ~ WIS d (x, a, 0, 1%; t), the formula above 
becomes even simpler and we have for v = vr + ivi such that vr G l\ a; t, vj G <Sd(R): 

E[exp W „X,-))] = •^y-yP. (") 
det(/ rf - 2tl^v) 2 
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1.3 Some identities in law for affine processes 

This section presents simple but interesting identities in law for affine processes. We first state a prelim- 
inary lemma. 

Lemma 6 — Let B G C(S d (R)) and q G d (R). We define B q G £(<S d (R)) 6y = (q T )- 1 B(q T xq)q- 1 . 

Then, 

B satisfies (|4]) <^=> £> g satisfies 

Proof : It is sufficient to prove the direct implication since B(x) = q T B q ((q T )^ 1 xq^ 1 )q. Let x, v G <5j~(R) 
such that Tr(a;w) = 0. We have 

Tr(B,(a» = Tr[ J B(g T x 9 ) 9 -^(? T )" 1 ] > 0, 
since B satisfies (HJ and r Vr[{q T )~ 1 xq~ x qvq T ] = Tr(xv) =0. □ 

Proposition 7 — The following identities holds: 

• AFFd(x,a,B,a) = AFF d (x, a, B, \ / a T a). 

Lau) 

• (linear transformation) Let q G £?d(R), we have 

q T AFF d (x,a,B,a)q = AFF d (q T xq,q T aq,B q -i,aq), 

Law 

where B q -i is defined by Vy £ 5d(R), B q -i(y) = q T B((q T )~ 1 yq~ 1 )q. 

Proof : From Proposition [3] or Corollary |H we know that (Xf) t >o ~ AFF d (x,a, B,a) and (X t 2: ) t >o ~ 
AFF d (x, a, B, V a 1 a) have the same infinitesimal generator. Therefore, they solve the same martingale 
problem. Thanks to the affine structure, one can show that this martingale problem has a unique solution 
by looking at the characteristic function. This is made in a much general case in Cuchiero et al. [5] . Similarly, 
we can check easily that q T AFF d (y, a, B, a)q and AFF d (q T yq, q T aq, B q -i, aq) lead to the same martingale 
problem. Here, we remark that B q -i satisfies ((H) thanks to Lemma [5] □ 

An interesting consequence of this linear transformation is given in the following corollary. It states that 
any affine process can be obtained as a linear transformation of an affine process for which we have a = 1%, 
where L^ = (^i=j<n)i<i,j<d (see notations in Appendix . Since our main goal here is to sample paths of 
such processes, this says us that it is sufficient to focus on this special case. 

Corollary 8 — Let (Xf) t >o ~ AFF d (x, a, B, a) and n = Rk(a) be the rank of a T a. Then, there exist a 
diagonal matrix 8, and a non singular matrix u G 0<j(R) such that a — u T 5u, and a T a = u T L^u, and we 
have: 

(A7) t > = u T AFF d {(u- 1 ) T xur' L ,8,B u ,I2)u, 

Law 

where \/y G S d (R), B u (y) = (u _1 ) T B(u T yu)u _1 . 
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Proof : Once u is given, the identity in law comes directly from Proposition [7] We give now a constructive 
proof of the existence of it, which takes back the arguments given by Golub and Van Loan (|12j. Theorem 
8.7.1). Nonetheless, we explain it entirely since it gives a practical way to get u. 

Let us consider a + a T a e 6>j~([R). From the extended Cholesky decomposition given in Lemma [33] there 
is a matrix v G Gd{^) such that v T av + v T a T av — L d , where r = Rk(a + a T a). Since v T av G <Sj~(R), 
v T a T av € <Sj~(IR) and z 1 I r d z = for z € R d such that z\ = ■ ■ ■ = z r = 0, there are si, S2 & $n (R) such that: 

T - fsi 0\ jTT / S 2 \ 

w av — I q I and v a av — I ^ o / ' 

Let 02 be an orthogonal matrix such that o^sioi is a diagonal matrix. We assume without loss of generality 
that only the first n elements of this diagonal are positive: 0^202 = diag(rjx, . . . , r) n , 0, . . . , 0). We set 

o = ( °q and get I r d — o T v T avo + o T v T a T avo, which gives that o T v T avo is a diagonal matrix. 

Thus, we get the desired result by taking u — diag(y/rji, . . . , y 7 ?^, 1, . . . , l)o~ 1 v~ 1 . □ 



Let us make few comments on the practical implementation to compute u. For Wishart processes, 
a — aa T a, and u can directly be obtained by using a single extended Cholesky decomposition (Lemma [33]) . 
In the general case, the computation of u mainly requires an extended Cholesky decomposition and the 
diagonalization of s 2 . 

Up to now, we have stated identities for the law of affine processes. Thanks to the explicit characteristic 
function of Wishart processes, we are also able to get another interesting identity on the marginal laws. 

Proposition 9 — Let t > 0, a, b € A^d(R) o,nd a > d— 1. Let m t — exp(ifr), q t — J Q exp(s&)a T a exp(sb T )ds 
and n — Rk(g t ). Then, there is 9t E Gd(R) such that q t — t0 t I d 0j ' , and we have: 

WIS d (x,a,b,a;t) = 9 t WIS d (9t 1 m t xm£(ft 1 ) T , a, 0, 7J; t)6? (12) 

Law 

This proposition plays a crucial role for the exact simulation of Wishart processes. Thanks to (fT2")l , we can 
sample any Wishart distribution if we are able to simulate exactly the distribution WLSd(x, a, 0, L d \ t) for 
any x € <Sj~(R). In Section[2] we focus on this and give a way to sample exactly WIS d (x, a, 0, L d ;t). Let us 
stress here that we can compute the matrix 0t by using the extended Cholesky decomposition of qt/t, as it 
is explained in the proof below. 

Proof : We apply Lemma l33l to qt/t € <S^ (R) and consider (p, c n , k n ) an extended Cholesky decomposi- 
tion of q t /t. We set 9 t = p^ 1 ( ? n T ^ ) . Then, Q t is invertible and it is easy to check that q t — tOtl^Oj '. 
Now, let us observe that for v S Sd(R), 

det(J d - 2iq t v) = dct(6» t (6' < r 1 - 2itl d l 6'[v)) = det(J d - 2itL d l 6fv6 t ), 
Tr[iv(I d - 2iq t v)~ 1 m t xmJ] = Tr[i(6» t " 1 ) T 6» t T u(6»t6i f r 1 - 2it6 t I26?v9 t 6t 1 )- 1 m t xmT} 

= Tr[i9?v9 t (I d - 2itI26jve t )~ 1 6- 1 m t xmJ{e- 1 ) T ]. 

Let X? ~ WIS d (x, a, b, a; t) and X? ~ WIS d (x, a, 0, /£; t). Then, from $TU§ and CCE]), we get that 

E[exp(rTr(vX?))] = E[exp(iTr(0f u0 t xf r 1 ™ m TV7 1 ) T ))] = E[exp(m(v9 t xt* lm ^V* 1)T 
which gives the result. □ 



(i 



Last, we have to mention that the identity (fT2|) extends an usual identity between CIR and Bessel squared 
distribution. It gives when d — 1: 

2bt _ i 2btx 

WISx(x,a,b,a\t) = a 2 —-—WIS 1 (-^ r = ^,a,0,l;t). 

Law lot a z (l — e zot ) 

In that case, this identity can also be obtained directly from the SDE. Let (X£ )t>o ~ WISi(x, a, b, a). We 
have dXf = (aa 2 + 2bX t )dt + 2a^/Xf dW t . Then, Y t = e~ 2bt Xf jo? is a time-changed Bessel squared process 

since dY t = a(e- 2bt dt)+2y/Y t (e- bt dW t ). We obtain WISi (x, a, b, a; t) = a 2 e 2bt WLSi{x/a 2 , a, 0, 1; )■ 

Law 

On the other and, a linear time-change gives that WISi(x, a, 0, 1; Xt) = XWIS\{x/ X, a, 0, 1; t), which leads 

Law 

to the same identity as {H} by taking A = (1 - e~ 2bt ) / (2bt) . 

2 Exact simulation of Wishart processes 

In this section, we present a new method to simulate exactly a Wishart process. It works without any 
restriction on the parameters. Wishart distributions have been thoroughly studied in statistics when a G N 
(which is then called the number of degrees of freedom). Exact simulation methods have already been 
proposed in that case. For instance, Odell and Feiveson [21] have given an exact simulation algorithm for 
central Wishart distributions based on the Bartlett's decomposition, and Gleser [TT| extends it to any (non- 
central) Wishart distribution. Bru [1] explains also when a € N how Wishart processes can be obtained as 
a square of Ornstein-Uhlenbeck processes on matrices. 

Here, our method relies on the identity in law (IT21 that enables us to focus on the case 6 = 0, a = FJ. 
Then, we show a remarkable splitting of the infinitesimal generator as the sum of commuting operators. 
These operators are associated to SDE that can be solved explicitly on iSj~([R), which enables us to sample 
any Wishart distribution. 

2.1 A remarkable splitting for WISd(x, a, 0, 1%) 

The following theorem explains how to split the infinitesimal generator of a WISd(x, a, 0, JJ) as the sum 
of commutative infinitesimal generators. It will play a crucial role in the sequel for the exact simulation and 
to get discretization schemes. 

Theorem 10 — Let L be the generator associated to the Wishart process WISd{x,a,0, 1%) and Li be the 
generator associated to WLSd{x, a, 0, e l d ) for i £ {1, . . . , d}. Then, we have 

n 

L = ~^^Liand Vi,j€.{l ) ...,d},LiLj=LjLi. (13) 

i=l 

Proof : From ©, we easily get that L = 53lLi since i? = Y^i=i e d- The commutativity property comes 
from a simple but tedious calculation which is left in Appendix IC.ll □ 

Beyond the commutativity, two other features of (|13|) are important to notice. 

• The operators L t and Lj are the same up to the exchange of coordinates i and j. 

• The processes WLSd(x,a,0 7 e d ) and WISd{x, a, 0, 7J) are well defined on <Sj~(IR) under the same hy- 
pothesis, namely a > d — 1 and x € Sj~(IR). 
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This second property enables us the composition that we explain now. Let us consider t > and x € <Sj~(IR). 
We define iteratively: 



WLS d (x,a,0,e d ;t), 



WIS d (Xl' x ,a,0,elt), 



Thus, X\'"' is sampled according to the distribution at time t of a Wishart process starting from 

i-i x t 1,a: 

X t '"' and with parameters (a, 0, e d ). We have the following result. 



Proposition 11 Let X™ 



be defined as above. Then 



WIS d (x,a,0,I2;t). 



Thanks to this proposition, we can generate a sample according to WIS d (x, a, 0, JJ; t) as soon as we can 
simulate the laws WIS d (x, a, 0, e d ). These laws are the same as WIS d (x, a, 0, I d ;t), up to the permutation 
of the first and i th coordinates. In the next subsection, it is explained how to draw such random variables. 

It is really easy to give a formal proof of Proposition [TTJ Let / be a smooth function on <Sj~(IR) and 
X? ~ WIS d (x, a, 0, 12; t). By iterating Ito's formula, we have that E[f{Xf)] = J2T=o t k L k f(x)/kl. Similarly, 
we also get by using the tower property of the conditional expectation that: 



E 



n, . . . 



= E 



E 



n-1,... t 



h n ■ 



E 

fe„=0 



L k nf{xr l,. x * ) 



(14) 



Simply by repeating this argument, we get that 

+ °° t E? =1 



E 



f(x n,... t ) 



; l L^...L k r rf(x)=J2j-,(Li + ---+L n ) k f(x) = E[f(Xf)]. (15) 
fci,...,fc„=0 1 71 ' fc=0 



E 



To get the second equality, we identify a Cauchy product and use that the operators L\, . . . ,L n commute. 
To make this formal proof correct, one has to check that the series are well defined and can be switched with 
the expectation. This check is made in the Appendix IC.2I for our framework and remains valid as soon as 
the operator Li and L are of affine type. 



2.2 Exact simulation for WIS d (x,a,0, I^;t) 

In this subsection, we focus on the simulation of WLS d (x, a, 0, /]) with a > d- 1 and x £ Sj(R). Thanks 
to PropositionfTTl this enable us then to simulate samples of WLS d (x, a, 0, L2). For the sake of the clearness, 
we start with the case of d = 2 that avoids complexities due to the matrix decompositions. We deal with 
the general case just after. 

2.2.1 The case d = 2 

We start by writing explicitly the infinitesimal generator Li of WIS2{x, a, 0, From ©, we get: 



x e S+(R), Ltf(x) = ad {1<iy f(x) + 2x {1 ^d\ lM f{x) + 2x {h2} d {1 . 1} d {1 . 2} f(x) + ^ld 2 {h2] f(x). (It 
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We show now that this operator is in fact associated to an SDE that can be explicitly solved. We will denote 
by {Z},t > 0) and (Zf,t > 0) two independent standard Brownian motions. 

When £{2.2} = 0, we also have £{1,2} = since x is nonnegative. In that case, 

Xg = x, ) {M} = adt + 2^{Xf) {lM dZl d{X? ) {1)2} = 0, d(Xf ) {2 , 2} - (17) 

has the infinitesimal generator (|16p . which is the one of a Cox-Ingersoll-Ross SDE (or of a squared Bessel 
process of dimension a to be more precise). By using an algorithm that samples exactly a non central 
chi-square distribution (see for instance Glasserman [10 ), we can then sample WIS%(x, a, 0, 1%; t) when 

£{2,2} = 0. 

When £{2,2} > 0, it easy to check that the SDE 



<W) { i,i> = a dt + 2^ (X ? ) {M} - { -^f^fdZj + * jgfe^l 

<*(*f){i,2> = v / (^f){2,2}^ 2 (18) 
TOW) = 



starting from Xq — x has an infinitesimal generator equal to L\. To solve (JTHJ) , we set: 

(TJU\ / vx\ 

i(Xt){l,2}) ( TT u\ (-^f ){1,2} / rr «\ 

(.^t ){i,i> = ( x t ){i,i} — (Y^r r — ~' \ u t hm - PF —— ' (^t ){2,2} - £{2,2}- 

l A t J{2,2} V X {2,2} 

Here, u stands for the initial condition, i.e. u = Uq. We get by using ltd calculus that 



(19) 



d(Un { iA} = (a - l)dt + 2^(U?) {lil} dZi, d(U?) {1 , 2} = dZ 2 t and d{U?) {m = 0. (20) 

Therefore, (f/"){i,2} and (£^"){i,i} can be respectively be sampled by independent Gaussian and non-central 
chi-square variables. Then, we can get back Xf by inverting (|19l) : 

(Xf) {1 , 1} = (U?) {1 , 1} + (f/ t ")f 1>2} , (*f){l,2} = (E?){l,2}^/(t?){2, 2 }, (Xn { 2,2} = {Ut){2.2}- (21) 

The following proposition sums up and makes more precise the above result. Its proof is postponed to 
the general case (Theorem [TJ 



Proposition 12 — Let x € S 2 (R). T/ien, i/ie process (Xf) f > defined by either (fTTf w/ien £{2,2} = 
or (|18[) w;/ien £{2,2} > ^ as ^ 5 infinitesimal generator equal to L\. Moreover, the SDE (|18[) /ias a unique 

strong solution which is given by (|21[) . where (U",t > 0) is £/ie solution of the SDE (|20[) starting from 

2 

"{1,1} - £{1,1} ~ sgff > 0, U {1)2} = ^==, ^{2,2} = £{2,2}- 

This result gives an interesting way to figure out the dynamic associated to the operator L\, by using 
a change of variable. It is interesting to notice that the CIR process (Ut){i,i} is well defined as soon 
as its degree a — 1 is nonnegative, which coincides with the condition under which the Wishart process 
WIS2(x, a, 0, J 2 ) is well-defined. Last, we notice that the solution of the operator L\ involves a CIR process 
in the diagonal term and a Brownian motion in the non diagonal one. We will get a similar structure in the 
general d case. 
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2.2.2 The general case when d > 2 

We present now a general way to sample exactly WISd(x,a,0,I^;t). We first write explicitly from © 
the infinitesimal generator of WISd{x, a, 0, 1]) for x € 5j"(R): 

^{i,m}i9{i,m}9{i,i}/(a;) + x X {m,l}d{l, m }d {1 J}f{x). 



l<m<o 



l<m,i<d 



(22) 

As for d = 2 we will construct an SDE that has the same infinitesimal generator L\ and that can be solved 
explicitly. To do so, we need however to use further matrix decomposition results. In the case d = 2, we 
have already noticed that we choose different SDEs whether £2,2 = or not. Here, the SDE will depend on 
the rank of the submatrix {xi.j)2<i,j<d, and we set: 

r = Rk((xi,j)2<i,j<d) G {0, . . . , d - I}. 
First, we consider the case where 



3c r £ Q r lower triangular, k r G .A/f £ ;_i_ rxr ([R), (x) 



2<i,j<d 



C r 
k r 







(23) 



With a slight abuse of notation, we consider that this decomposition also holds when r = with c = 0. 
When r = d — 1, c = c r is simply the usual Cholesky decomposition of {xij)2<i,j<d- As it is explained in 
Corollary Q21 we can still get such a decomposition up to a permutation of the coordinates {2, . . . , d}. 



Theorem 13 — Let us consider x € 5j~(R) such that (|23j) holds. Let [Z\~)i<i< T +i be a vector of independent 
standard Brownian motions. Then, the following SDE (convention Efc=i(- ■ • ) = when r = 0) 



d(Xf) {hl} = adt + 2 V /(*f) {1 , 1} - ELi (E[ =1 (cr X )m(*? ){!,!+!}) d^t 1 

+2ELi E[=i(c7 1 ) fc ,K^f){i,/ +1 }< +1 (24) 



<*(*f){i,<} 

^((^f){i,fe})2<fc,i<d 



Efc=i Ci-i,kdZ t 




fe+i 



2,...,d 



/ias a unique strong solution (Af) t >o starting from x. It takes values in 5j~([R) and has the infinitesimal 
generator L\. Moreover, this solution is given explicitly by: 



X? 



1 

C r 

k r I(i— r —\ 



( 



(U?){1,1} + E ((W){l,fc+1}) ((un{i,i + i})l<i<r 
fe=l 

((Ut) {1,1 + 1}) l<l<r It 

/ 



1 

c r k r 

I d -r-l 



(25) 



whe 



d (( C/ "){l,/+l})l<i<r 



(a - r)dt + 2 v /(C/ t ") {ljl} dZ t 1 , = ^{1,1} - ELi( u {i,fc+i}) 2 ^ °' (26) 

(dZ f )i<;< r , (u{ij+i})i<;<r = c ^ 1 (a ; {i,;+i})i<;<r- 



Af = 



Once again, we have made a slight abuse of notation when r = 0, and (|25p should be simply read as 
I in that case. In the statement above, it may seem weird that we use for u and 











the same indexation as the one for symmetric matrices while we only use its first row (or column). The 
reason is that we can in fact see Xf as a function of Z7" by setting: 

i U t){i,j} = u {i,j} = x {i,j} for 3 > 2 and (E/"){i,i} = m {m} = for r + 1 < i < d. (27) 
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Thus, (c r , k r , Id-i) is an extended Cholesky decomposition of ((U™)ij)2<ij<d and can be seen as a function 
of Ul 1 . We get from (gSJ that 



h r {u) 



1 

Cr(u) 
k r {u) Id-r-l 



X? = h(U?), with h{u) = Y^%=vk[{u i , j h< itS < d ]hr{u) and 

/ r \ 

«{i,i} + E ( w {i,fe+i}) 2 )i<;<r 

V 



(28) 



1 

c r (u) T k r {u) T 

I d - r -i 



k=l 

( u {i,;+i})i<;<r ir 

/ 

where (c r (u), k r {u), Id-i) is the extended Cholesky decomposition of (lfj,j)2<i,j<d given by some algorithm 
(e.g. Golub and Van Loan [TJ], Algorithm 4.2.4). Equation (|28p will play later an important role to analyse 
discretization schemes. 

The proof of Theorem [TjJ] is given in Appendix IC.31 It enables us to simulate exactly the distribution 
WISd(x, a, 0, I\; t) simply by sampling one non-central chi-square distribution for {Ul')^^ (see Glasser- 
man [10j ) and r other independent Gaussian random variables. Like in the d — 2 case, we notice that 
the condition which ensures that the Cox-Ingersoll-Ross process ((L r "){i,i} , i > 0) is well defined for any 
r € {0, . . . , d—1}, namely a— (d— 1) > 0, is the same as the one required for the definition oiWISd(x, a, 0, 



Remark 14 — From (|25p . we get easily by a calculation made in (|47[) that Y(k(Xf) — Rk((xij)2<ij<d) + 
1(C/«) {1 1} ^o, and therefore, 

Rk(Xf) = Rk((x lJ ) 2 <. i , J < (i ) + 1, a.s. 

Theorem [T3l assumes that the initial value a; <G 5j"(R) satisfies (|23|) . Now, we explain why it is still 
possible up to a permutation of the coordinates to be in such a case. This relies on the extended Cholesky 
decomposition which is stated in Lemma l33l 



Corollary 15 — Let (Af) t >o ~ WISd(x, a, 0, l\) and (c r ,k r ,p) be an extended Cholesky decomposition of 
(xi,j)2<i.j<d (Lemma \33\) . Then, 7r= fg ^ J is a permutation matrix, and 



{Xt)t>a = TT T WIS d (nxiT T ,a, 0,/])tt and ((■Kxir T ) itj ) 2 < i . <d = ( 

law \ 



C r 
k r 







satisfies 



Proof : The result is a direct implication from Proposition [71 since 7r T = 7r 1 and ~kI\-k t = l\ 



□ 



Therefore, by a combination of Corollary [15] and Theorem ll3[ we get a simple way to construct explicitly 
a process that has the infinitesimal generator L\ for any initial condition x £ S^(R). In particular, this 
enables us to sample exactly the Wishart distribution WISd{x, a, 0, l\; t). The algorithm below sums up 
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the whole procedure. 



Algorithm 1: Exact simulation for the operator L\ 



Set 7T 



Input : x e 5j"(R), d, a > d - 1 and t > 0. 

Output: X, sampled according to WISd(x, a, 0, l\;t) 

Compute the extended Cholesky decomposition (p,k r ,c r ) of (xi,j)2<i,j<d given by Lemma 133] 
r e {0, . . . , d — 1} (see Golub and Van Loan [12] for an algorithm); 

q p \ x = nxir T , (u{i,/+i})i<i< r = (cr) _1 (*{i,i+i})i<i<r and 

= ^{1,1} - ELi( w {i,fc+i}) 2 ^ ; 

Sample independently r normal variables G2, . . . , G r+ i ~ 7V(0, 1) and (£/"){!,!} as a CIR process at 
time i starting from u.n n solving d([/ t u ){i.i} = (a — r)eft + 2^(U^)^i A jdZ^ (See Glasserman [TO])- 
Set (f/ t u ) {lji+1} = + VtGi+i ; 

return X = 




/ 



- E ((t^){i,fc+i}) 

t{l,i + l})l<!<r 





((^ t u ) { 



l,i+l}Jl<i<'- 









/ 











Id-r-1 



Let us discuss now the complexity of this algorithm. The number of operations required by the extended 
Cholesky decomposition is of order 0(d 3 ). From a computational point of view, the permutation is handled 
directly and does not require any matrix multiplication so that we can consider w.l.o.g. that ir = Id- Since 
c r is lower triangular, the calculation of i = 1, ...,r + 1 only requires 0{d 2 ) operations. Also, we 

do not perform in practice the matrix product (1251) . but only compute the values of for i = 1, . . . ,d, 

which requires also 0(d 2 ) operations. Last, d samples are at most required. To sum up, it comes out that 
the complexity of the whole algorithm is of order 0(d 3 ). 

2.3 Exact simulation for Wishart processes 

We have now shown all the mathematical results that enable us to give an exact simulation method 
for general Wishart processes. This is made in two steps. First, by using the remarkable splitting (Propo- 
sition [IT]) and the exact scheme for WISd{x, a, 0, I^;t) (Theorem IT51 and Corollary [T5]) . we get an exact 
simulation scheme for WISd(x, a, 0, 12; t). As we will see, it is related but extends the Bartlett's decompo- 
sition of central Wishart distribution that dates back to 1933. Then, by using the identity in law (IT21) . we 
are able to sample any Wishart distribution WISd{x 1 a, b, a; t). 

2.3.1 Exact simulation for WIS d (x, a, 0, 7J; t) 

By gathering the results obtained in the subsections 12.11 and 12.21 we get a way to sample exactly the 
distribution WISd(x,a,0,I2;t). Indeed, thanks to Theorem [TBI and Corollary [TBI we know how to sample 
exactly WISd(x, a, 0, 1^; t). By a simple permutation of the first and k th coordinates, we are then also able 
to sample according to WISd(x,a,0,e^;t) for k € {!,..., d}. Thus, we get by Proposition ITTI an exact 
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simulation method to sample WISd(x, a, 0, J?; t). It is given explicitly in the algorithm below. 

Algorithm 2: Exact simulation for WISd{x, a, 0, JJ; t) 

Input: x € 5j~(R), n < d, a > d — 1 and £ > 0. 
Output: X, sampled according to WISd{x, a, 0, I^; t) 
y = x 

for fc = 1 to n do 

Set pk,i — Pi.k — Pi,i = 1 for i $ {1 5 k}, and pj j = otherwise (permutation of the first and k th 
coordinates). 

y = pYp where Y is sampled according to WISd{pyp, a, 0, I\; t) by using Algorithm [T] 
end 

return X = y. 

Since this algorithm basically runs n times Algorithm [TJ it requires a complexity of order 0(nd 3 ) and 
therefore at most of order 0(d A ). As we have seen, the "bottleneck" of Algorithm Q] is the extended Cholesky 
decomposition which is in 0(d 3 ). All the other steps in Algorithm [1] require at most 0(d 2 ) operations. A 
natural question for Algorithm [2] is to wonder if we can reuse the Cholesky decomposition between the for 
loops instead of calculating it from scratch. For example, if it were possible to get the Cholesky decomposition 
of loop k+ 1 from the one of loop A; at a cost 0(d 2 ), the complexity of Algorithm [5] would then drop to 
0(d 3 ). Despite our investigations, we have not been able to do so up to now. 

Remark 16 — When a > 2d — 1, it is possible to sample WISd{x, a, 0, I^', t) in 0(d 3 ) by another mean. If 
X\ - WIS d (x, d, 0, I^; t) and X? - WIS d (0, a - d, 0, 1%; t) are independent, we can check that X\ + Xf ~ 
WISd{x, a, 0, J?; t). Then, X\ can be sampled by using Provosition \30\ and X 2 by using Corollary \18\ since 
Xl = tWIS d (0,a-d,0,I2;l) from $n$. 

Law 



Remark 17 — Let x £ Sj~ : *([R). From Remark^TQ we get that for each k e {!,..., n}, X t '"' * € 
iSj~'*(IR), a.s.. In that case, the extended Cholesky decomposition that has to be computed at each step is an 
usual Cholesky decomposition. 



2.3.2 The Bartlett's decomposition revisited 

Now, we would like to illustrate our exact simulation method on the particular case WISd(0, a, 0, 1%; 1), 
which is known in the literature as the central Wishart distribution. In that case, we can perform explicitly 
the composition given by Proposition [TT] We will show by an induction on n that: 

v n,... x i" _ ( (Li,j)l<i,j< n \ / {Ljj)l<i,j< n \ 

Xl - ) \ ) ' 

where {L.^^Kj^^d and L^j are independent random variables such that Li j ~ A/"(0, 1) and (L^.i) 2 ~ 
X 2 (oi — i + and L t j — for i < j. 

For n = 1, we know from Theorem [TBI that (Xl'°)i t i ~ X 2 (a) since d(Xl'°) 1A = adt + 2^{Xl fl ) lA dZl 

with {X Q ' )x i = 0, and all the other elements are equal to 0. Let us assume now that the induction hypothesis 
is satisfied for n — 1. Then, we can apply once again Theorem 1131 (up to the permutation of the first and 

x i,o 

n th coordinates). We have Rk(X"~ '"' ) = n — l,a.s., and the Cholesky decomposition is directly given 
by (Li,j)i<i,j<n-i- Then, we get from (|25[) that there are independent variables L 2 n ~ X 2 ( a — ti + 1) and 
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1,0 

n,../'l 



L n ,i ~ Af(0, 1) for i G {1, . . . ,n - 1} such that X[ 

l 

10)1 (-£<n,<)l<j<n-l E2=l^,< I I 10 

o o / d _„ / V o o o/V o o i d -, 











In-1 


(Ln,i)l<i<n-1 


1 







.)T 

,VKi<ri-l 


V" T 2 





Id-n 












(L n ,i) 


\<i<n-l 






1 




T 2 


:)■( 


(-^n : i)i<i<„_i 











~ 



In-1 \I<n,i)l<i<n— 1 \ / -?n-l \ / J n _i (L„ ,)i<,<n-l 







(^i;j)j<i,j- 


:) 


( 
























) ( t 











Ln,n 



we conclude by induction and get the following result which is known as the Bartlett's decomposition (see 
Kshirsagar [17] or Kabe [IB]). 



Corollary 18 Bartlett's Decomposition of central Wishart distributions. With the notations 
above, 

( ;/ - ! ;, r 1 S)( ( ^ }l - w - n ° Q )~wis d (o,a,o,i2;i)- 

Therefore, the identity given by Proposition [II] that we have obtained by a remarkable splitting of 
the infinitesimal generator extends in a natural way the Bartlett's decomposition to non-central Wishart 
distributions. 



2.3.3 Exact simulation for WISd{x,a,b,a;t) 

Now, thanks to the identity in law (|12[) . we get an exact simulation scheme for WISd(x, a, &, a; t). The 
associated algorithm is given below. To the best of our knowledge, this is the first exact simulation method 
for non-central Wishart distributions that works for any a > d — 1. The existing methods in the literature 
mainly focus on the case a € N that arises naturally in statistics. Namely, Odell and Feiveson [21] and 
Smith and Hocking [35] have proposed an exact simulation method for central Wishart distributions based 
on the Bartlett's decomposition. Gleser [IT] has given an exact simulation method for non-central Wishart 
distributions, also when the degree a is an integer. 

Algorithm 3: Exact simulation for WISd(x, a, b, a; t) 

Input: x G Sj{R), a > d - 1, a, b G M d (R) and t > 0. 
Output: X, sampled according to WISd(x, a, b, a; t). 

Calculate qt = J Q exp(s6)a T a exp(s6 T )ds and (p, c„, k n ) an extended Cholesky decomposition of qt/t. 
Set 9t = p^ 1 ( f™ T ^ ) and m t = exp(tb). 

\ Kn ld—n J 

return X = t Y6j ', where Y - WISd{0^ 1 m t xmf(6- 1 ) T , a, 0, /J; t) is sampled by Algorithm^ 



3 High order discretization schemes for Wishart and semi definite 
positive affine processes 

Up to know we have given an exact simulation method for Wishart processes. It relies on a remarkable 
splitting of the infinitesimal generator given in Theorem 1101 When dealing with discretization schemes, 
splitting operators is a powerful technique to construct schemes for SDEs from other schemes obtained 
on simpler SDEs. This idea of splitting originates from the seminal work of Strang [23] in the field of 
ODEs. As pointed by Ninomiya and Victoir [20] or Alfonsi [I], it rather easy to analyse the weak error 
of schemes obtained by splitting, simply by using the same arguments as Talay and Tubaro [24) for the 
Euler-Maruyama scheme. Thus, we will only focus on the analysis of the weak error, i.e. the error made on 
marginal distributions. 
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Up to our knowledge, there are very few papers in the literature that deal with discretization schemes 
for Wishart processes. Recently, Benabid et al. [2] have proposed a Monte-Carlo method to calculate expec- 
tations on Wishart processes which is based on a Girsanov change of probability. Gauthier and Possamai [3] 
introduce a moment-matching scheme for Wishart processes. Both methods are well defined under some re- 
strictions on the parameters, and there is no theoretical result on their accuracy. Currently, Teichmann |25j 
is working on dedicated schemes for general affine processes by approximating their characteristic functions. 

This section is structured as follows. First, we recall basic results on the splitting technique to get 
discretization schemes for SDEs. We will take the same framework as Alfonsi PQ since it is somehow designed 
for affine processes. Then, we will explain how to get high order schemes for WISd{x, a, 0, 7J) from the 
remarkable splitting (TT5]) . From this result, we will be able to get a second order scheme for any semi 
definite positive affine processes and a third order scheme for Wishart processes. 

3.1 General results on discretization schemes for SDEs 

In this paragraph, we try to present the minimal knowledge that we will use later to get schemes for 
affine processes. We take back the framework developed in Alfonsi ([T], Section 1) and refer to this paper 
for further details and proofs. 

Let us start with some notations. We consider a domain D C R\ ( £ N*. For a multi-index 7 = 
(71, ... ,7{) € we define 9 7 = dj 1 , . . . , d^ c and I7I = Yli=i 7» an d set: 

Cprf(ID) = {/ G C°°(D, R), V 7 G N c , 3C 7 > 0, e 7 G IN*, Vx G D, \d y f(x)\ < C 7 (l + ||x|| e ^)}, 

where ||.|| is a norm on R^. We will say that (C 7 ,e 7 ) 76N c is a good sequence for / G C^j(D) if one has 
|<9 7 /(x)| < C 7 (l + ||cc|| e7 ). Mainly (but not only), we will consider in this paper D = <Sj"(R) as a subspace 
of S d (R) ~ R d ( d+1 )/ 2 . In that case, it is natural to keep the same indexation as for matrices, and we rather 
use the notation 8 1 = Ux<i<j<d d {lj} for 7 = (7{i,i})i<i<i<d G N d ( d+1 )/ 2 . 

Definition 19 — Let b : D -> R c , a : D -> M((R). The operator L defined for f G C 2 (D, R) by: 

Lf{x) =J2 b i^)9if(x) + - (aa T ) ilj (x)d i d j f(x) (29) 

is said to satisfy the required assumptions on if the following conditions hold: 
. Vi, j G {1, . . . , C}, bi(x), (o-a T ) id {x) G C- ol (D), 

• for any x € D, the SDE Xf — x + J* b(X*)ds + f Q cr(Xf )dW s has a unique weak solution defined for 
t > 0, and thus P(Vt > 0, X'f € D) = 1. 

In the case of affine diffusions, bi(x) and (aa T )i.j (x) are affine functions of x and the operator satisfies the 
required assumption on the appropriated domain. Let us stress that if L satisfies the required assumption 
on D and / G C™ ol (D), all the iterated functions L k f(x) are well defined on D and belong to C^ ol (D) for any 
k G N. 

Let us turn now to discretization schemes. We will consider in this paper a final time horizon T > and 
the regular time grid defined by tf — iT/N, i = 0, . . . , N. When discretizing a Markovian process, a scheme 
is usually described as a way to sample the process at the next time step from its current position. Thus, we 
will say that (p x (t)(dz),t > 0,x G D) is a family of transition probabilities on D if p x (t) is a probability law 
on D for any time-step t > and any position x G D. Then, a discretization scheme on the regular time grid 
associated to this family starting from xo G D is simply a sequence {X]% , < i < N) of D-valued random 
variables such that: 

• X^ N = xq, 
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• the law of Xj^, is sampled according to p<y N (T/N)(dz) independently from the previous samples, i.e. 
E[f(X!i )\(X% ,0<j<i)]= L f(z)PxN (T/N)(dz) for any bounded measurable function / : D -> R. 

i+1 j t N 

Up to the initial condition, the discretization scheme (Xj%, , < i < N) is entirely characterized by (p x (t)(dz),t > 

0, x £ D). With a slight abuse of language, we will then say that (p x (t)(dz), t > 0, x £ D) is a scheme on D. 
We will denote by Xf a random variable on D which is sampled according to p x (t) . 



Definition 20 — Let L be an operator that satisfies the required assumption on D. A scheme {p x {t){dz) 1 1 > 
0, x £ D) is a potential weak vth-order scheme for the operator L if for any function f £ C^ ol (E)) with a good 
sequence (C 7 , e 7 ) 7 g N < , there exist positive constants C, E, and rj depending only on (C 7 ,e 7 ) 7eN < such that 



Vie (0,7?), 



k=l 



fix) 



<Cr +1 (l + ||a:|r). 



(30) 



When the coefficients of the SDE have a sublinear growth, we can check that the exact scheme (i.e. Xf = Xf , 
where (Xf,t > 0) solves the SDE associated to L) is a potential ^th-order scheme for any order v £ tKl ([IJ, 
Proposition 1.12). In that case, (|30[) is then clearly equivalent to: 



3C,E,r)>0, Vt£ (0, 77), E[f(Xf)]-E[f(Xf)} < Ct^(l + \\x\\» ). 

In practice, having a potential weak ^th-order scheme for the operator L is the main requirement to get 
a weak error of order v. This is precised by the following theorem that relies on the idea developed by Talay 
and Tubaro l24l. 



Theorem 21 — Let us consider an operator L that satisfies the required assumptions on D and a discretiza- 
tion scheme (X/^,,0 < i < N) with transition probabilities p x (t)(dz) on D that starts from Xj% = xq £ D. 
We assume that 

1- Px{t){dz) is a potential weak vth-order scheme for the operator L, 
2. the scheme has bounded moments, i.e.: 

Vg £ N*, 3iV(q) £ N, sup E[||l t ^|| 9 ] < oo, (31) 

N>N(q),0<i<N i 



3. f : D — >• R is a function such that u(t,x) — E[f(X^_ t )] is defined on [0,T] x D, C°° , solves Vt £ 
[0,T],Va; £ D, dtu(t, x) — —Lu(t,x), and satisfies: 

VI £ N, 7 £ N c ,aC*/, 7 ,e/, 7 > 0,Vxe D,t € [0,T], \d l t d jU (t,x)\ < C*/, 7 (l + ||x|| e ^). (32) 
Then, there is K > 0, N £ N, such that \E[f(X$)] - E[f(X*°)]\ < K/N v for N > N . 

Let us mention that condition (13"TT) is slightly weakened with respect to Theorem 1.9 in [T]. However, we can 
check that ([3Tj) is in fact sufficient to make work the proof of this theorem given in [1] . 

Let us comment briefly the hypothesis of this theorem. The boundedness of the moments (f3"Tj) holds as 
soon as the coefficients of the SDE have a sublinear growth (see Lemma [351 in Appendix [Dj . It is instead 
much more technical to get the bounds (1321) on the derivatives of u. In their original paper, Talay and 
Tubaro have considered stronger assumptions on b and a (namely, C°° with bounded derivatives) to get such 
bounds. In the case of Wishart processes, we are able to get (|32"]) when / £ C^ ol (Sd{R)) by using the explicit 
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formula of the characteristic function (|TU| . This is stated in Proposition H31 We however deem that this 
result still holds for the general affine case. 

Now, we present results that explain how the property of being a potential weak ^th-order scheme can 
be preserved when we use the splitting technique. To do so, we need to introduce first the composition of 
schemes. Let us consider two transition probabilities p\{t){dz) and p 2 ,(t)(dz) on D. Then, we define: 

p 2 {t 2 )opl{t 1 ){dz) := f pKhXdz^ihXdy), 

which is the law obtained when one uses first use scheme 1 with a time step t\ and then scheme 2 with a time 
step t2 with independent samples. More generally, if one has m transition probabilities p\, . . . , on D, we 
define 

P m (t m ) o • • ■opl(t 1 )(dz) := P m (t m ) o (p m -\t m -l) o • • ■ apl(ti)(dz)). 



Proposition 22 — Let Li,L 2 be two operators satisfying the required assumption on D. Let p\ and p\ he 
respectively two potential weak vth-order schemes on D for L\ and L 2 - 

• If L1L2 — LiL\, p 2 (t) o p\{t){dz) is a potential weak vth-order discretization scheme for L\ + h%. 

• If v > 2, p 2 (t/2) °j5 1 (t) op^(t/2) and i (p 2 (t) ° p\{t) + p 1 (t) p 2 (t)) are potential weak second order 
schemes for L\ + L 2 - 



3.2 High order schemes for Wishart processes 

In this paragraph, we will give a way to get weak i/th-order schemes for any Wishart processes. The 
construction of these schemes is the same as the one used for the exact scheme. First we obtain a i/th-order 
scheme for WISd(x, a, 0,23). Then, we get a z/th-order scheme for WISd(x, a, 0, I d l ) by scheme composition. 
Last, we use the identity in law (fT2j) to get a weak z/th-order scheme scheme for any Wishart processes. 



3.2.1 A potential weak i/th-order scheme associated to L\ 

In this paragraph, we present a potential weak ^th-order scheme for the generator of WISd(x, a, 0,1^). 
Roughly speaking, we obtain this scheme from the exact scheme given by Theorem [13] and Corollary [15] by 
replacing the Gaussian random variables with moment matching variables and the exact CIR distribution 
with a sample according to a potential weak v-th order scheme for the CIR. 



Theorem 23 — Let x G <Sj~(R) and (c r ,k r ,p) be an extended Cholesky decomposition of (%i,j)2<i,j<d- We 

set 7T = ( \ ^ I and x = irxir T , so that x — ( f r !? I I C I ^ I ■ Like in Theorem \1SI we have 
\ p J \ k r \ ) \ 10 \ ) 

r 

"{1,1} = £{M} - ^2( u {i,k+i}) 2 > 0, where (u { i j;+ i})i<;< r = c^ 1 (x {1 .i+i})i<i< r , 
fe=i 

and we set it{i,i} = Oi/r + 2 < i < d and uuj\ = xuj\ if ' i,j > 2. Let (G l )i<i< r be a sequence of 
independent real variables with finite moments of any order such that: 

Vi€{l,...,r}, Vfc< 2u + l, E[{G l ) k ] = E[G k ], where G ~ W(0, 1). 

Let (Ut){i,i} be sampled independently according to a potential weak vth-order scheme for the CIR process 
d^U^^nj — (a — r)dt + 2^J (U^^i^jdZ^ starting from M{i.i} and that satisfies the immersion property 
(Definition^. We set: 



(Pt){i,i} = HiM + ^\ 2 < i < r + 1, (l> t u ) {M } = 0, r + 2 < i < d, (U?) {itJ} = u {i , j} ifi,j > 2. 
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Then, the scheme defined by Xf = tt t h r (U^)TT is a potential vth-order scheme for L\, where the function 
h r is defined by (|28[) . 

The proof of this result is technical and requires further definitions. It is left in Appendix 1231 Second and 
third order schemes for the CIR process that satisfy the immersion property can be found in Alfonsi (T] (see 
Corollary |40|) . We can therefore get second (resp. third) order schemes for L\ by taking any variables that 
matches the five (resp. the seven) first moments of Af(0, 1). This can be obtained by taking 

P((7 = ^3) = P(G l = -VS) = - and P(G l = 0) = - (33) 

6 3 

(resp. p(&=ey/3 + V6) = P (& = e^3 - VG \= \ - ee{-l,l}). (34) 

3.2.2 A potential weak i^th-order scheme associated to WIS<i(x, a, 0, JJ) 

From TheoremfTUl we know that the infinitesimal generator of WISd(x, a, 0, J?) is given by L = J^ILr 
where the operators Li are the same as L\ up to the permutation of the first and i th coordinate. Moreover, 
these operators are commuting. Thus, if ir 1 ^ 1 denotes the associated permutation matrix and Xf is the 
potential i/th-order scheme defined by Theorem [23} 

^l-o-tj^-Tr 1 **^ ++ * 7r i++« j g a p t en ti a l i/th-order scheme for Lj. 
We then get a potential j/th-order scheme for L thanks to Proposition |2"2"1 

Corollary 24 — Let p l x , i G {1, . . . ,d} be potential weak vth-order scheme for Li. Then for any n < d, 

p n (t)o...opl(t)(dz) 

defines a potential weak vth-order scheme for WISd{x, a, 0, 1%). 

Remark 25 — Let e > and Xf be a potential weak vth-order scheme for WISd(x,a,0, Since 
x-\-et v+l Id is a potential weak vth-order scheme for the operator L = 0, we easily get by scheme composition 
(Proposition^^ that X^ +et Id is also a potential weak vth-order scheme for WISd(x, a, 0, 1%). This 
scheme starts from an invertible initial condition, and by Remark \14\ we only make in that case usual 
Cholesky decompositions on invertible matrices. 

3.2.3 A weak ^th-order scheme for Wishart processes WISd(x, a, b, a) 

Now that we have a potential ;/th-order scheme for WISd(x, a, 0, 7? ), we are in position to construct 
a scheme for any Wishart process WISd(x, a, b, a) thanks to the identity (fl2|). Unfortunately, we need to 
make some technical restrictions on a and b (namely, a 6 Gdffi) or ba T a — a T ab) to show that we get like 
this a potential ^th-order scheme. We however believe that this is rather due to our analysis of the error 
and that the scheme converges as well without this restriction. We mention in addition that we give in the 
next section a second order scheme based on Corollary [8] for which we can make our error analysis for any 
parameters. 

Proposition 26 — Lett > 0, a, b e A / Jd(R) and a > d—1. Let nit = exp(tb), q t = f Q exp(sb)a T aexp(sb T )ds 
and n — Rk(a T a). We assume that either a £ Gdffi) or b and a T a commute. We define 

• if n — d, Q t O.S the (usual) Cholesky decomposition ofqt/t, 
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• if n < d, 9 t = J j J* exp(sfe) exp(sb T )dsp 1 ( f" ^ J where (c n ,k n ,p) is the extended Cholesky 



decomposition of a a otherwise. 

In both cases, 9 t G Gd{R)- Let Yj denote a potential weak vth-order scheme for WISd(y,a,0, /?). Then, 
the scheme defined by 

X? = e t Yy mtXm * (e * 1)T 9?, (35) 
is a potential weak vth-order scheme for WISd(x, a, &, a). 

Proof : First, let us check that 9 t e Gdffi) is we U dehned, such that q t /t = 9 t I d l 9j an d satisfies: 

3K,r) >0, Vte (O.^.maxdl^lUl^ir 1 ) < K. (36) 

When n = d, qt/t is definite positive as a convex combination of definite positive matrices and the usual 
Cholesky decomposition is well defined. Moreover, (|36p holds since qt/t goes to a T a which is invertible 
when t — > + . When n < d, we have assumed in addition that b and a T a commute. Therefore, qt = 
a T a(J Q exp(sb) exp(sb T )ds/t). Since a T a and (J exp(sb) exp(sb T )ds/t) are positive semidefinite matrices 
that commute, we have 



qt = y — J exp(sb)exp(sb T )dsa T a^J - J exp(sb) exp(sb T )ds. 
Once again, j J* exp(s6) exp(sb T )ds is definite positive as a convex combination of definite positive matri- 



c„ 



ces and we get that 9 t = J \ J exp(sb) exp(sb T )dsp 1 f , I G Gd{R) satisfies q t /t = 9 t I d l 9f by 

Lemmal33l Similarly. (|36|) holds since p~ 1 ( f" ^ ) does not depend on t and i/ \ f* exp(s6) exp(s6 T )ds 

\ tin ld—n J " 

goes to Id when M0 + . 

Let / £ C^ oi (S/[ (R)). Let Xf ~ WISd{x, a, 6, a; t). Since the exact scheme is a potential i/th-order 
scheme, there are constants E,r/ > depending only on a good sequence of / such that 

Vt G (0,77), - £ ^/(z)l < ^ +1 (1 + INI*)- (37) 

fc=0 

On the other hand we have from Proposition [SJ 

E[/(Xf)] - £[f(Xf)] = E[f(6 t Yy mtXm t {e7l)T 9j)} - E[f{e t Y t 9 * lmtxm * {e * 1)T 9f)]. (38) 

Let us introduce fe t (y) := f(8 t y9f) G C™ ol (S£(R)). By the chain rule, we have %j}/0 t (y) = Tr[6» t (e^ + 

^]e J /)9fdf(8 t y9f)}, where {df(x)) k ,i = (1a=i + ^O^fe,*}/^) and e 1 / = (l fc=i , l=3 -)i<M<d- From ® , 
we see that there is a good sequence (C 7 , e 7 ) 7gN d(d+i)/2 that can be obtained from a good sequence of / such 
that: 

Vie (0,r/),VyG«S+(R), |a y /* t (y)|<C T (l + ||y|| e -'). 
Therefore, we get that there are constants still denoted by C, E, r\ > such that 

vt g (0,7?), E[f(9 t Y t e " mtxmT{e " )T 9f)} - E[/(0 t y/ rWro * T(e <~ 1)T t T )] < +1 (i + ii^Wmn^fll 5 )- 

(39) 

From we get that there is a constant K' > such that ||6' t " 1 m t a;rn^(6' t " 1 )' r || £; < if' ||a;|| B for t € (0,7/). 
Thus, we get the result by gathering J37|, (J3HJ) and (JS3) - □ 



19 



Now, we want to show thanks to Theorem [5T] that the scheme given by Proposition [55] gives indeed a 
weak error of order v. However, if we exclude the exact simulation of the CIR, we only have at our disposal 
(up to our knowledge) at most a third order scheme for the CIR (see pQ). It seems thus fair to state the 
result on the weak error in that case. 

Theorem 27 — Let (Xf) t >o ~ WIS d (x, a, b, a) such that either a £ Gd(R) or a T ab — ba T a, and f £ 
C™ ol (S d (R)). Let (X$,0< i< N) be sampled with the sch erne defined by Provosition \26\ with the third order 
scheme for the CIR given in JJ^ and starting from x £ <Sj~([R). Then, 

3C,N > 0,VN > N a , \E[f(X t %)} - E[f(X?)}\ < C/N 3 . 
Proof : The conditions of Theorem [5T] are satisfied thanks to Propositions [551 EH an d Lemma [351 a 



3.3 Second order schemes for affine diffusions on <Sj~(IR) 

In this part, we present two different potential second order schemes for general affine processes. The 
first one is well-defined without any restriction on the parameters. The second one is faster but requires to 
assume in addition that a — da T a £ S^(R). 

Let (Xf)t>o ~ AFFd(x, a, B, a). Thanks to CorollaryO there is u £ Gd(R) and a diagonal matrix S such 
that a — u T Su, a 7 a = u T I^u and we have: 

(X x )t>o = (u T Y t {u " )Txu "u) t > , where Y t y ~ AFF d (y,S,B u , 1%). 

law 



Lemma 28 — If Y t y is a potential vth-order scheme for AFFd(y,6, B u , 1^), then u T Y^ u ^ xu u is a 
potential vth-order scheme for AFFd(x,a, B, a). 

Proof : Let / £ C™ ol (Sj(R)). We then have x i— > f(u T xu) £ C™ ol (S+(R)). Since u is fixed, there are 
constants C,7],E depending only on a good sequence of / such that for t £ (0,1]), \£[f(u T Y} u ' xu u)] — 

E[f(xf)}\ = |E[/( u T f/ u_1)Tra_1 u )] - w\f{u T Y^ 1)Txu ^u)}\ < cr +1 {i + iKu- 1 )^- 1 ]^) < ct» +1 {i + 

\\x\\ E ), for some constant C > C. □ 

We now focus on finding a scheme for AFFd(y,a,B u ,I2). Since 6 is a diagonal matrix such that 
8-(d- 1)12 G S d( R )> we have 

5 min := min 5^ >d-l. 

l<i<n 

From Corollary [4l we can write the infinitesimal generator of Y t v 
L = Tr([£ + B(x)]D s ) + 2Tv(xD s I2D s ) = Tv([S - S min L d l + B u {x)]D s ) + 5 m - m Tr{D s ) + 2Tv(xD s ' I^D 8 ) 

V v ' V v ' 

Lode I J wis d (x,s Jain ,o,i2) 

(40) 

as the sum of the infinitesimal generator of WISd{x, <5 m i n , 0, 1%) and of the generator of the affine ODE 
x'it) = S — 5minI2 + B u (x(t)). Of course, this ODE can be solved explicitly. Moreover, we know by 
Lemma l45l that if x(0) — x £ St(R), the solution of this ODE remains in 5j~(IR) since Assumption (|4| holds 
for B u and 5 — <5 m i n -T^ £ Sf(R). We denote by p^ DE (t) the Dirac mass at x(t): this is an exact scheme for 
Lode- Thanks to Proposition [22j we get the following result. 
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Proposition 29 — Let p^(t) denote a potential second order scheme for WISd(x, S m i n , 0, 1^). With the 
notations above, the following schemes 

\p° DE {t)°€{t) + \p W (t)op° DE (t) andp ODE (t/2)op^(t)opO^ {t /2) 

are potential second order scheme for the affine process AFF d (x, S, B u , I^). 

By combining Lcmma[28]and Proposition[29l we get a potential second order scheme for any affine process. 
In the numerical experiments in Section IU we have used the composition p ODE {t/2) o p w (t) o p® DE (t/2) 
even though the other one would have work as well. 

Let us stress now that different splitting from (|4TJ)) are possible. For example, we could have chosen instead 
L = Tr([5-pi2 + B u (x)}D s )+l3Tr(D s ) + 2Tr(xD s I2D s ) for any /3 G [d-l,S min ]. When a-da T a G <S+(R) 
(or equivalently when S — dl% G 6>j~(R)), the following splitting 

L = Tr([£ - dI2 + B u (x)]D s ) + dTr(D s ) + 2Tr(xD s F l fD s ) (41) 

V V 

Lode ^wiSj^.d.ci") 

is really interesting. Indeed it is known from Bru [I] that Wishart processes can be seen as the square of an 
Ornstein-Uhlenbeck process on matrices and can be simulated very efficiently More precisely, we will use 
the following result. 

Proposition 30 — Let x G <Sj~(R) and c G Md(R) be such that c T c — x. Then, the process Xf = 
(c + WJ2) T (c + W f L d l ) satisfies 

(X?)t> = WIS d (x,d,0J d l ). 

law 

If G denote a d-by-d matrix with independent elements sampled according to (|33j> . Xf = (c+ ^JtGI d l ) T {c + 
\ftGI~2) is a potential second order scheme for WISd(x, d, 0, JJ). 

Proof : We have by using Ito calculus dXf = (c + W t l d ) T dW t l d + I2dW?{c + W t I$) + dl^dt. By 
using Lemma The quadratic covariation of (Xf)ij and [Xf) m ^ n is given by d{{Xf (Xf ) m , n ) — 
(X?)i im (I2) jt n + {X*) i>n {I2) j<m + (XfhmWhn + ( X t km- Therefore, {Xf) t > solves the'same 

martingale problem as WISd(x, d, 0, 12), which is known to have a unique solution from Cuchiero et al. [BJ. 

Let us show now that Xf is a potential second order scheme. We can see c + \ftGL2 as the Ninomiya- 
Victoir scheme with moment- matching variables (see Theorem 1.18 in [1) associated to | Xa=i Xw=i 
on M d (R). Let / G C™ ol (S% (R)) . Then, x G M d (R) i-> f(x T x) G C™ ol (M d (R)) and there are constants 
C, E, ij > depending only on a good sequence of / such that: 

V* G (0, r,), |E[/((c + VtGI2) T (c + VtGI2))} - E[/((c + W t l2) T (c + W t l2))}\ < Ct» +1 (l + \\c\\ E ). 

Let us observe now that the Frobenius norm of c is ^Tr(c T c) = v/Tr^J < y/d + Tr(x 2 ) < Vd + ^Ty(x 2 ). 
Therefore, for any norm, there is a constant K > such that ||c|| < K(l + \\x\\), which gives the result. □ 



Corollary 31 — We assume that 6 — dI2 G <Sj~(R). Let p Wbi3 (t) x denote the second order scheme of 
Provosition \3(A and p® DE (t) be the exact scheme for the ODE associated to Tr([<5 — dI2 + B u (x)]D s ). Then, 
\p ODE (t) opf»-{t) + \p w »"{t) o p° DE (t) andp ODE {t/2) o p w »-(t) o p° DE {t j2) are potential second-order 
scheme for AFF d (x,8,B u ,I2). 
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The scheme given by Lemma |2"51 and Corollary [3"T1 is not well-defined when a — da T a £ <5>j~([R). When it is 
defined, it is however more efficient than the scheme given by Proposition[29]and Lcmma[28] We have already 
mentioned in subsection 12.3.11 that the exact scheme has a cost of 0(d 4 ). Similarly, the schemes given by 
Propostions [2"o1 and [21)1 have a cost of 0(d 4 ) operations. On the contrary, the scheme given by Propositionl3"01 
only requires one Cholesky decomposition. Thus, the scheme given by Lemma [2H1 and Corollary [3T] has a 
complexity in 0(d 3 ). We will illustrate this in the next section. 

We conclude this section by mentioning that in the case of Wishart processes, we can show that the weak 
error is indeed of order 2 thanks to Theorem |2"T1 Lemma 131)1 and Proposition H31 

Theorem 32 — Let {Xf) t > ~ WIS d (x,a,b,a). Let (X^,0< i < N) be the sch erne starting from x which 
is sampled either by the scheme defined by Lemma \28\ and Proposition \29\ ( with the second order scheme for 
the CIR given in f^) or, if we have in addition a > d, by the scheme given by Lemma\E^ and Corollarv \31\ 
Then, 

3C,N Q > 0,ViV > N , \E[f(X^)] - E[/(Xf)]| < C/N 2 . 



4 Numerical results on the simulation methods 

The scope of this section is to compare the different simulation methods given in this paper. We still 
consider a time horizon T and the regular time-grid if = iT/N , for i — 0, . . . , N. In addition, we want 
to compare our schemes to a standard one, and we will consider the following corrected Euler-Maruyama 
scheme for AFFd(x,a,B,a): 

X$ = x, X$ = X*n + (a+B(X$ ) ) — + y/(X£ )+(W t N +i -W t N)a+a T (W t N +i -W t N) T ^(X^)+,0 < i < N-l. 

(42) 

Here, x + denotes the matrix that has the same eigenvectors as x with the same eigenvalue if it is positive 
and a zero eigenvalue otherwise. Namely, we set x + — odiag(\f , . . . , A^)o T for x — odiag{\\, . . . , \d)o T . 
Thus, x + is by construction a positive semidefinite matrix and its square root is well defined. Without this 
positive part, the scheme above is not well defined for any realization of W . 

First, we compare the time required by the different schemes and the exact simulation. Then, we present 
numerical results on the convergence of the different schemes. Last, we give an application of our scheme to 
the Gourieroux-Sufana model in finance. 

4.1 Time comparison between the different algorithms 

In this paragraph, we compare the time required by the different schemes given in this paper. As it has 
already been mentioned, the complexity of the exact scheme as well as the one of the second order scheme 
(given by Lemma [2"51 and Proposition [2^]) and the third order scheme (given by Proposition [25]) is in 0(rf 4 ) 
for one time-step. To be more precise, they require 0{d ) operations that mainly corresponds to d Cholesky 
decompositions, 0(d 2 ) generations of Gaussian (or moment-matching) variables and 0(d) generations of 
nonccntral chi-square distributions (or second or third order schemes for the CIR). The time saved by the 
second and third order schemes with respect to the exact scheme only comes from the generation of random 
variables. For example, the generation of the moment-matching variables ([33"]) and (f34|) is 2.5 faster than 
the generation of A/"(0, 1) on our computer. The gain between the second or third order schemes for the 
CIR given in Alfonsi [1] and the exact sampling of the CIR given by Glasserman [10 is much greater, but it 
depends on the parameters of the CIR. When the dimension d gets larger, the absolute gain in time between 
the discretization schemes and the exact scheme is of course increased. However, the relative gain instead 
decreases to 1, because more and more time is devoted to matrix operations and Cholesky decompositions 
that are the same in both cases. Let us now quickly analyse the complexity of the other schemes. The 
second order scheme given by Lemma [28] and Corollary [31] (called "second order bis" later) has a complexity 
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N = 10 


N = 30 


Schemes 


R. value 


Im. value 


Time 


R. value 


Im. value 


Time 


Exact (1 step) 


-0.526852 


-0.227962 


12 








2 nd order bis 


-0.526229 


-0.228663 


41 


-0.526486 


-0.229078 


125 


2 nd order 


-0.526577 


-0.228923 


76 


-0.526574 


-0.228133 


229 


3 rd order 


-0.527021 


-0.227286 


82 


-0.527613 


-0.228376 


244 


Exact (N steps) 


-0.526963 


-0.228303 


123 


-0.526891 


-0.227729 


369 


Corrected Euler 


-0.525627* 


-0.233863* 


225 


-0.525638* 


-0.231449* 


687 


a = 3.5, d = 3, A B , = 10" 


3 ,A /m = 10- 


~ 3 , exact value R. = — 


3.527090 and Im.= -0.228251 


Exact (1 step) 


-0.591579 


-0.037651 


12 








2 nd Qrder 


-0.590444 


-0.037024 


77 


-0.590808 


-0.036487 


229 


3 rd order 


-0.591234 


-0.034847 


82 


-0.590818 


-0.036210 


246 


Exact [N steps) 


-0.591169 


-0.036618 


174 


-0.592145 


-0.037411 


920 


Corrected Euler 


-0.589735* 


-0.042002* 


223 


-0.590079* 


-0.039937* 


680 


a =2.2, d = 3, A R = 0.9 x 10" 


3 , Aj m = 1.3 x 10 3 , exact value R. 


= -0.591411 and Im.= - 


0.036346 


Exact (1 step) 


0.062712 


-0.063757 


181 








2 nd order bis 


0.064237 


-0.063825 


921 


0.064573 


-0.062747 


2762 


2 nd Qrder 


0.064922 


-0.064103 


1431 


0.063534 


-0.063280 


4283 


3 rd order 


0.064620 


-0.064543 


1446 


0.064120 


-0.063122 


4343 


Exact (N steps) 


0.063418 


-0.064636 


1806 


0.063469 


-0.064380 


5408 


Corrected Euler 


0.068298* 


-0.058491* 


2312 


0.061732* 


-0.056882* 


7113 


a = 10.5, d = 10, A R = 1.4 x 10" 3 , A Im = 1.3 x 10" 3 , exact value R. = 0.063960 and Im.= - 


0.063544 


Exact (1 step) 


-0.036869 


-0.094156 


177 








2 nd Qrder 


-0.036246 


-0.094196 


1430 


-0.035944 


-0.092770 


4285 


3 rd order 


-0.035408 


-0.093479 


1441 


-0.036277 


-0.093178 


4327 


Exact (N steps) 


-0.036478 


-0.092860 


1866 


-0.036145 


-0.093003 


6385 


Corrected Euler 


-0.028685* 


-0.094281* 


2321 


-0.030118* 


-0.088988* 


7144 



a = 9.2, d = 10, A R = 1.4 x 10 3 , A Im = 1.4 x 10~ 3 , exact value R. = -0.036064 and Im.= -0.093275 



Table 1: E[exp(— TriivX^r))] calculated by a Monte- Carlo with 10 6 samples for a Wishart process with a = Id, 
b — 0, x — 10Id, v = 0.09/d and T = 1. The starred numbers are those for which the exact value is outside the 95% 
confidence interval, and Ar (resp. Ai) gives the two standard deviations value on the real (resp. imaginary) part. 

in 0(d 3 ) operations for one Cholesky decomposition and matrix multiplications, with 0(d 2 ) generations of 
Gaussian variables. The complexity of the corrected Euler scheme is of the same kind. At each time-step, 
0(d 3 ) operations are needed for matrix multiplications and for diagonalizing the matrix in order to compute 
the square-root of its positive part. However, diagonalizing a symmetric matrix is in practice much longer 
than computing a Cholesky decomposition even though both algorithms are in 0(d 3 ). Also, one has to 
sample 0(d 2 ) Gaussian variables for the Brownian increments. 

In Table 14.11 we have calculated by a Monte-Carlo method one value of the characteristic function of a 
Wishart process. It is also known analytically thanks to (|10l) . and we have indicated in each case the exact 
value. We have considered dimensions d — 3 and d — 10. We have given in each case an example where 
a > d and another one where d — 1 < a < d. We have used the different algorithms presented in this paper: 
"2 nd order bis" stands for the scheme given by Lemma [25] and Corollary [21] (with the moment-matching 
variables ([3"3"]) h "2 nd order" stands for the scheme given by Lemma |2"51 and Proposition [23] (with ([3"3"| and the 
second order scheme for the CIR given by pQ), "3 rd order" stands for the scheme given by Proposition [26] 
(with ([34]) and the third order scheme for the CIR given by [TJ), and "Corrected Euler" stands for the 
corrected Euler-Maruyama scheme presented in the beginning of this section. For the exact scheme, we have 
both considered the cases with one time-step T and N time-steps T/N. Of course, the first case is sufficient 
to calculate an expectation that only depends on Xt, but the second case allows to compute also pathwisc 
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expectations. For each method, we have given the value obtained and the time needed on our computer 
(3000 MHz CPU). 

First, let us mention that the exact value is in each case in the confidence interval except for the corrected 
Euler scheme. As one can expect, the exact method with one time-step is by far the quickest method to 
compute an expectation that only depends on the final value. We put aside this case and focus now on the 
generation of the whole path. We see from Table 14.11 that the second and the third order schemes require 
roughly the same computation time. As expected, the second order scheme bis is much faster when it is 
defined (i.e. when a> d). On the contrary, the Euler scheme is much slower than the second and third order 
scheme, even though it should be faster for large d. This is due to the cost of the matrix diagonalization. 
Let us mention that the time required by the discretization schemes is proportional to N and do not depend 
on the parameters when the dimension is given. On the contrary, the time needed by the exact scheme 
may change according to a and can increase considerably when a is close to d — 1. To be more precise, the 
exact simulation method for the CIR given by Glasserman [10] uses a rejection sampling when the degree 
of freedom is lower than 1, which corresponds to the case d — 1 < a < d. The rejection rate can in fact be 
rather high, notably when the time-step gets smaller. For N = 30, d = 3 and a — 2.2, the exact scheme is 
four times slower than the second order scheme and 2.5 slower than the exact scheme with a = 3.5. 

Let us draw a conclusion from this time comparison between the different schemes. Obviously, we 
recommend to use the exact scheme when calculating expectations that depend on one or few dates. Instead, 
when calculating pathwise expectations of affine processes by Monte-Carlo, we would recommend to use in 
general the second order bis scheme when a > d and the second order (or third order for Wishart processes) 
when d — 1 < a < d. 

4.2 Numerical results on the convergence 

Now, we want to illustrate the theoretical results of convergence obtained in this paper for the different 
schemes. To do so, we have plotted for each scheme E[exp(—Tr(ivX^))] in function of the time step T/N. 
This expectation is calculated by a Monte-Carlo method. As for the time comparison, we illustrate the 
convergence for d = 3 in Figure Q] and d = 10 in Figure O Each time, we consider a case where a > d and a 
case where d — 1 < a < d, which is in general tougher. In these figures, 

• scheme 1 denotes the value obtained by the exact scheme with one time-step, 

• scheme 2 stands for the second order scheme given by Lemma [2FJ and Proposition [251 

• scheme 3 denotes the third order scheme given by Theorem |2"TI 

• scheme 4 is the corrected Euler scheme (|4"2"j) . 

Here, we have not plotted the convergence of the second order (bis) scheme given by Lemma [25] and Corol- 
lary [21] because it would have given almost the same convergence as the other second order scheme. 



N 


2 


4 


8 


10 


16 


30 


Figure [T] right 
Figure [2] right 


-0.000698 
0.494752 


0.000394 
-0.464121 


0.033193 
0.657041 


0.111991 
0.643042 


0.185128 
0.637585 


0.210201 
0.619553 



Table 2: Values obtained by the Euler scheme in the numerical experiments of Figures Q] and [2] 

As expected, we observe in both Figures Q] and [2] convergences that fit our theoretical results. Namely, 
Scheme 2 converges in 0(1/ N 2 ) and Scheme 3 converges faster in 0(l/iV 3 ). In some cases such as Figure [2] 
Scheme 3 already matches the exact value from N = 2. Even though it seems to converge at a 0(1/N) 
speed, the corrected Euler scheme is clearly not competitive with respect to the other schemes. In the tough 
case d — 1 < a < d, the values obtained by the Euler scheme are in fact outside the figures, and we have put 
the corresponding values in Tabled] 
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Figure 1: d = 3, 10 7 Monte-Carlo samples, T = 10. The real value of E[exp(— Tr(ivX^,))] in function of the time-step 
T/N. Left: v = 0.05/d and Wishart parameters x = 0AI d , a = 4.5, a = I d and 6 = 0. Exact value: 0.054277. Right: 
v = 0.2I d + 0.04q and Wishart parameters x = 0.4/ d +0.2q, a = 2.22, a = I d and b = ~0.5I d . Exact value: 0.239836. 
Here, q is the matrix defined by: qi.j — li^j. The width of each point represents the 95% confidence interval. 

We want to conclude this section by testing numerically the convergence of our schemes when we cal- 
culate pathwise expectations. Of course, our theoretical results only bring on the weak error, but we 
may hope that our schemes converge also quickly when considering more intricate expectations. In Fig- 
ure [4?2j we approximate E[maxo<t<T Tr(Xf)] with the different schemes by computing the maximum on 
the time-grid. The convergence seems to be roughly in 0(1 /V~N) for all the schemes (see Figure B~2l 
left), including the exact scheme. However, the main error seems to come from the approximation of 
maxo<t<r TifXT ) by maxo<fc<2V Tr(Xwv). In fact, we have plotted in Figure l4~2l (right) the difference be- 

— — — — t k 

tween E[maxo<fc<iV Tr(Xi^)] and E[maxo<fc<jv Tr^ 2 ^)]. Then, we find convergences that are very similar 
to those obtained for the weak error: schemes 2 and 3 converge at a speed which is respectively compatible 
with 0(1/7V 2 ) and 0(1/N 3 ). Scheme 4 seems also to give a 0(1/N) convergence. It would be hasty to 
draw a global conclusion from this simple example. Nonetheless, the convergence of schemes 2 and 3 is 
really encouraging on pathwise expectations, if we put aside the problem of approximating a function of 
{X? , < t < T) by a function of (Xf N , < k < N). 

4.3 An application in finance to the Gourieroux and Sufana model 

In this paragraph, we want to give a possible application ot our schemes in finance. More precisely, 
we will consider the model introduced by Gourieroux and Sufana 13J. This is a model for d risky assets 
Si , . . . , Sf . Let (B t , t > 0) denote a standard Brownian motion on R d that is independent from (Wt,t > 0). 
Then, we consider the following dynamics for the assets: 

t > 0, 1 < I < d, S l f = S l + r f S l u du + f S l u {^X~ u dB u ) u (43) 

Jo Jo 

where X t = X a + J* (aa T a + bX u + X u b T ) du + f* (y^X^dW u a + a T dW^^/X^) is a Wishart process. Here, 
(V X u dB u )i is simply the I th coordinates of the vector \/X u dB u . We can easily check that the instantaneous 
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1 2 3 4 5 6 7 8 9 10 



Figure 2: d = 10, 10 7 Monte-Carlo samples, T = 10. Left: imaginary value of E[exp(-Tr(ivX^ N ))] with v = 0.009/ d 
in function of the time-step T/N. Wishart parameters: x = OAId, a — 12.5, b — and a = Id- Exact value: 

-0.361586. Right: real value of Efexpf-TrfwX^))] with v = 0.009/ d in function of T/N. Wishart parameters: 

N 

x — OAId, a — 9.2, b = —0.57^ and a = Id- Exact value 0.572241. The width of each point represents the 95% 
confidence interval. 



quadratic covariation matrix between the log-prices of the assets is X t . Last, r denotes the instantaneous 
interest rate. 

To simulate both assets and the Wishart matrix, we proceed as follows. We observe that the generator 
of (St,X t ) can be written as 

d 1 ^ 

L = L s + L x , where L s = ^ rsid St + - ^ s,s r r, j). i). , 

i—l — \ 

and L x is the generator of the Wishart process WISd{x, a, b, a). The operator L s is associated to the SDE 
dS\ = rS\ + S l t (^/xdB t )i that can be solved explicitly. We have indeed S l t = Sq exp[(r — xij/2)t + {\fxB t )i\. 
Let us also remark that \fxB t = cB t if we have cc T — x: both are centred Gaussian vectors with the 

Law 

same covariance matrix. In practice, it is more efficient to use S\ = S^exp^r — xij/2)t + (cB t )i] where 
c is computed with an extended Cholesky decomposition of x rather than calculating y/x, which requires 
a diagonalization. Let us then denote by pj? x Jt) the exact scheme for L° (x is unchanged) and x Jt) 
a second order scheme for the Wishart process WISd(x,a,b,a) (s is unchanged). Then, we consider the 
following scheme 

\{p X (t) °Pfs, x) (t) +P S (t)°P( SlX )(t)), 

i.e. we draw a Bernoulli variable of parameter 1/2 to decide whether we use first p^ s x ^(t) or pj? x ^(t). Under 
some conditions that we do not check here, this construction is known from Proposition [22] to preserve the 
second-order convergence. To be consistent with Subsection 14.21 this scheme is denoted by scheme 2 later 
in this paragraph. To compare this scheme with a more basic one, we consider the Euler-Maruyama scheme 
defined by (J42]) and 

S[f = S l , ^jf = S l t f (l + tT/N + {JxJ f {B t N +i ~ B t? )),) ,0 < i < N - 1. 
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Figure 3: d = 3, 10 7 Monte-Carlo samples, T — 1. Wishart parameters x = OAId + 0.2g with q t j — li^j, a — 2.2, 
6 = and a — Id- Left, E[maxo<fc<iv Tr(X^f)], right: E[max <fc<jv Tr(X^,)] — E[max <fc<jv Tr(X* N )] in function of 
T/N. The width of each point gives the precision up to two standard deviations. 

It is denoted by scheme 4 like in Subsection 14.21 

We have plotted in Figure l4~3l the price of a put option on the maximum of two risky assets (d — 2). The 
Gourieroux and Sufana model is an affine model, and the characteristic function of St is explicitly known 
(see [13]). Thus, it is possible to adapt the method proposed by Carr and Madan [5] and to calculate by 
numerical integration (which is possible for small dimensions) the value of this put option. We have given 
in Figure 14.31 the exact value obtained by this method. As one might have guessed, we observe a quadratic 
convergence for scheme 2 and a linear convergence for scheme 4. The benefit of using scheme 2 is clear since 
it already fits with the exact value from N — 5 in both cases: its convergence is really satisfactory. 

5 Conclusion and prospects 

Let us draw a brief summary of this paper. Thanks to a remarkable splitting of the infinitesimal generator 
of Wishart processes, we have been able to sample exactly any Wishart distribution. We have also proposed 
a third order scheme for Wishart processes and a second order scheme for general affine diffusions. We have 
confirmed these rates of convergence with numerical tests and analysed the time complexity of each method. 
It comes out that we recommend to use the exact scheme to compute expectations that depend on one 
(or few) times. To calculate pathwise expectations, we instead recommend generally to use discretization 
schemes. More precisely, the second order scheme given by Lemma [2"51 and Corollary [31] has to be preferred 
when a > d. Otherwise, we recommend to use the third order scheme given by Theorem [23 for Wishart 
processes or the second order scheme given by Lemma [28] and Proposition 1291 for general affine diffusions. 

Let us give now some prospects of this work. As a possible continuation of this paper, it is natural to 
wonder if it is possible to extend our schemes to affine diffusions on positive semidefinite matrices that include 
jumps (see Cuchiero et al.[6]). From a modelling point of view, we believe that Wishart processes could be 
used in a wide range of applications. In fact, they can be used as soon as one has to model dependence 
dynamics. Thus, we hope that the possibility of sampling such processes will stimulate different kinds of 
dependence models. 
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Figure 4: E[e~ rT (K -maxiS 1 '* , S 2 N N )) + ] in function olT/N. d = 2, T = 1, K = 120, So = Sg = 100, and r = 0.02. 
Wishart parameters: x = 0.047 d + 0.02g with = 1^, a = 0.2I d , b = 0.5/d and a = 4.5 (left), a = 1.05 (right). 
The width of each point gives the precision up to two standard deviations (10 6 Monte-Carlo samples). 
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A Notations and some results on matrices 



Notations for real matrices : 

• For d, d' £ N*, .Md(R) denotes the real d square matrices and Mdxd'ffi) the real matrices with d rows 
and d! columns. 

• iSrf(IR), <Sj~(R), <Sj~'*(R), and ^(R) denote respectively the set of symmetric, symmetric positive semidef- 
inite, symmetric positive definite and non singular matrices. 

• For x £ Aid(R), x T , adj(x), det(x), Tr(x) and Rk(a;) are respectively the transpose, the adjugate, the 
determinant, the trace and the rank of x. 

• For x £ <Sj~(R), \fx denotes the unique symmetric positive semidefinite matrix such that (y/x) 2 = x 

• The identity matrix is denoted by Id and we set for n < d, 1% = {li=j< n )i<i,j<d and eJJ = (li=j =n )i<i,j<d, 
so that 12 = Ya=i e d- We also set for 1 <h3 <d, e% d = (^k=i,l=j)i<k,l<d- 

• For x £ 5d(R), we denote by the value of Xij, so that x = J2i<i<j<d x {i,j} ( e d 3 + : "-^i e d' i )- We use 
both notations in the paper: notation (:Eij)i<tj<d is of course more convenient for matrix calculations 
while (x{ijy)i<i<j<d is preferred to emphasize that we work on symmetric matrices and that we have 

X^j — Xj^i- 

• For Ai, . . . , Xd £ R, diag(X\, . . . , Xd) denotes the diagonal matrix such that diag(Xi, . . . , Xd)i.i = Xi. 

Now, we present different results that are used in the paper. We first recall the extended Cholcsky 
decomposition of positive semidefinite matrices which is used intensively and then give two other technical 
results. 

Lemma 33 — Let q £ iSj~ (R) be a matrix with rank r. Then there is a permutation matrix p, an invertible 
lower triangular matrix c r £ Gr(R) and k r £ ■Md-rxr(R) such that: 

T T ( C r \ 

MP = CC ' C= [kr J' 

The triplet (c r , k r ,p) is called an extended Cholesky decomposition ofq. Besides, c=(f r T ^ ) £ Gd(R), 
and we have: 

q = (c T p) T I r d C T p. 

The proof of this result and a numerical procedure to get such a decomposition can be found in Golub 
and Van Loan ([12], Algorithm 4.2.4). When r — d, we can take p = Id, and c r is the usual Cholcsky 
decomposition. 

Lemma 34 — Let b,c £ <S>d(R). If either b £ iSj~(R) or c £ S d (R), then Id + ibc is invertible. In particular, 
if b £ S d '*(R), b + ic is invertible. 

Proof : We start with the first assertion. Since (Id + ibc) T = Id + icb, it is sufficient to check the case 
where c £ Sj~(R). By a way of contradiction, let us assume that there is x £ C d \ {0} such that x + ibex = 0. 
We respectively denote by xr £ R d and xi £ R d the real and imaginary part of x. One gets easily that 
xji = bcxi and xj — —bcxR. Since x ^ 0, we have necessarily xr ^ 0, cxr ^ 0, bexn ^ and cbcxR ^ 0. 
Since c is nonnegative, we get by decomposing on an orthonormal basis that cxr.xr > and cbcxR.bcxR > 0. 
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However, we also have cxr.xr — —cxR.bcbcxR, which leads to a contradiction. The second assertion is now 
obvious since b + ic= b(I d + ib~ 1 c). □ 



B Proofs of Section [T] 

B.l Proof of Lemma [2] 

Proof : From the SDE ([5]) , we have 

(dY t ) id = (C t ) itj dt + (B t dW t A t + AfdW^Bf)^ 

= (Ct)ijdt + £fc )i=1 (CBtkfcUtkj + (Bt)jMAt)u)(dW t )kA 

Then we get the following formula for the quadratic covariation 

d{(Y t ) id , (Yt) m ,„) = Yf k ,i{{Bt)i,k{A t \j + {Bt) ,k{At)iMBt) m ,k{At)i,n + {B t ) n , k {At)i, m )dt 

= J2l,k(-Bt)i,k(Bt)m,k(A t )l, n {At)lJ + J2l,k(Bt)i,k(Bt)n,k(A t )l ym (A t )Lj 

+ Ef fe ( B t) J ,fc(SO m ,fc(A t ) i ,„(A t ) M +Ef, fc (-Btkfc(^)n,fc(^)i^(At) i , J (44) 
= {B t B t )i t7n {Aj A t ) n ,j + {BtBj)i^ n {Aj ' A t ) m ^ 

+ (B t B r [)j :m (Aj A t ) n ,i + (BtBf)j yn (Aj A t )m : i. 



B.2 Proof of Proposition [5] 

Proof : Let v £ 6>d(R) such that Vs £ [0, i\,I<i — 2q s v £ <?<z(IR). As it is usual for afhne diffusions, the 
Laplace transform can be formulated with ODE solutions. Namely, we will show that E[exp(Tr(vXf ))] = 
exp[4>(t, v) + Tr(/0(i, v)x)], where ip and solve following ODEs (see for example Cuchiero et al. [B]): 

d t tp(t,v) = ip(t,v)b + b T ip(t,v) + 2ip(t,v)a T cnp(t,v) ; ^(0,«) = v 
8t<t>it,v) = aTr(V>(*,u)) ; 4>(0, v) = 0. 

The function ip solves an usual matrix Riccati ODE. As shown by Levin |18j . ip can be obtained explicitly 
by the mean of the following exponential matrix : 

/ A M (t) A h2 (t) \ _ ( ( b -2a T a\\_( exp{tb) -2 exp(tb) f* exp(-sb)a T a exp(-sb T )ds) 
\A 2 ,i(t) A 3 , a (t) J 6XP \ \ -b T ))-{ exp(-i& T ) 

From Q2],we obtain that 

$(t,v) = (TP(0,v)A lt2 (t)+A 2 , 2 (t))- 1 (iP(0,v)A 1A {t)+A 2A (t)) 

= (^—2vexp(tb)J^exp(—sb)a T aexp(—sb T )ds + exp(—tb T )j vexp(tb) 

= (^—2vJ^exp(sb)a T aexp(sb T )dsexp(—tb T ) + cxp(—tb T )j vcxp(tb) 
= exp(tb T )(I d -2vq t )~ 1 vexp(tb), 

provided that I d — 2q s v is invertible for s £ [0, t], which holds by assumption. Therefore we get for x £ 5d(R), 

Tr(i/}(t,v)x) = Ti(exp{tb T ){I d ~2vqt)- 1 vexp(tb)x) 
= Tr ((I d - 2vq t )- 1 vexp(tb)xexp(tb T )) 
= Tr (v(I d - 2q t v)- 1 cxp(tb)xexp{tb T )) , 
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since v(Id — 2q t v) 1 = (Id ~ 2vq t ) 1 v. As explained by Grasselli and Tebaldi ([13], section 4.2), 4> can also 
be calculated explicitly by the mean of the exponential matrix above, and we get: 

0(t, v) = -|Tr (log[(7 d - 2w ft ) exp(i6 T )] - *&(&)) . 

By using that exp(Tr(log(A))) = det(A) for A G 0d(R), we deduce then that 

exp(0(t, w)) = exp(§ m(b)) (exp(Tr(log{(7 d - 2vq t ) exp(t6 T )})J) ^ 
= exp(f f&(6)) (det{(7 d - 2vq t )}det{exp{tb T )})^ 
= exp - f f&(6)) (det{(7 d - 2«<? t )})^ 

det(J d -2<j t i;)f ' 

For the last equality, we have used that det(ic/ — 2vqt) = det((id — 2vqt) T ) = det(7d — 2q t v). 
Now, it remains to show that (fT0|) indeed holds. By ltd calculus, we get that for s G (0,t): 

dexp[0(£-s, v)+Tr(ip(t-s, v)Xf)} = exp[0(t-s, u)+Tr(V>(i-s, u)X s x )]Trh/>(i-s, v)(v^dW s a+a T dW'J-v/*F)]- 

(45) 

Thus, exp[(j)(t — s,v) + Tr(^(t — s, v)Xg)] is a positive local martingale and therefore a supermartingale, 
which gives that E[exp(Tr(t>Xf ))] < exp[<p(t,v) + Tr(ip(t, v)x)] < oo, i.e. V b , a . t C 'D x ,a,b,a;t, where 

V b .a, t := {v G 5 d (R),Vs G [Q,t],I d - 2q s v G £ d (R)} and T> x , a ,b,a;t := G S d (R), E[exp(Tr(«X t x ))] < oo}. 

On the other hand, when — v G Sj~'*(R), we can check that exp[(j)(t — s,v) + Tr(ip(t — s,v)X^)] < 1 by 
observing that det(id — 2q t v) — det(7<2 + 2^/^vq t \/^v) > 1 and Tr (y(Id — 2q t v)~ 1 exp(tb)x exp(tb T )) = 
-Tr (^/^v(I d + 2y/^q t V^)~ 1 V^vexp(tb)xexp(tb T )) < 0. In that case, exp[<j>(t-s,v)+Tr('tp(t-s,v)Xf)] 
is a martingale from (|45|) . and (flOl) holds. 

Let us observe now that T>b,a-,t is convex. In fact, we have det(7d — 2q s v) = det(7d — 2 y /q^vy / q^), and 
therefore V bja . t = {d£ S^R), Vs G [0,t], 7<j — 2^/q^v^/q^ G Sj~'*(R)} which is obviously convex. The Laplace 
transform v H> E[exp(Tr(t?Xf ))] is an analytic function on 2\ a; t (see for example Lemma 10.8 in [8]). 
The RHS of (fT0|) is also analytic on T>b^ a -t and coincides with the Laplace transform when — v G <Sj~'*(R). 
Therefore, (ITU1) holds for v G T>b,a;t since 2\ a ;t is convex. Now, we can extend to complex values of v. Indeed, 
the RHS of (flOf is well defined for v = vr + ivi with vr G T> b ,a-,ti thanks to Lemma l34l Since both hand 
sides are analytic functions of v, (|10j) holds for v — Vr + iv%. 

Last, we want to show that T> b ^ a . t = Dx,a.b,a-.t- We first consider the case b = and assume by a way of 
contradiction that there is v G 'D x ,a,o,a;t \ T^o,a-t for some x, a, a and t > 0. Let i = min{s G [0, i] , I d — 2q s v ^ 
<?d(R)} G (0,t], On the one hand, we have v T> a .j and v G T> , a . s for s G [0, t). On the other hand, we 
have by Jensen's inequality: 

sG [0,i], exp((t- s)Tr(va T a))exp(Tr(vXf)) < E[exp(Tr(t)Xf ))|J r s ], 

which gives s G [Q,i] i-» exp(— sTr(ua T a))E[exp(Tr(wXf ))] is nondecreasing and finite. Since <jTU]) holds for 
s < t, we get that E[exp(Tr(wX?))] = +oo, which leads to a contradiction. Let us now consider the case 
b 0. From Proposition^ (which is a consequence of the characteristic function obtained above), we have 

V G V x , a ,b,a;t ^ 0?v9 t G V ,i S ;t ^ Vs G [0,t],det(7 d - 2{s/t)q t v) ? 0. 

In particular, f> x , a ,b,a;t is an open set. For v G £/d(R), we have det(7d — 2(s/t)q t v) ^ <^=> dct(t) -1 — 
2(s/t)q t ) 7^ (resp. det(7<j — 2g s u) ^ <^=> det(w~ 1 — 2q s ) ^ 0). Since sgt < s'qt (resp. g s < 
q s r) for s < s', we know from Theorem 8.1.5 in [12] that the (real) eigenvalues of v^ 1 — 2(s/t)qt (resp. 
v^ 1 — 2q s ) are nonincreasing w.r.t. s. Since they are also continuous, and v^ 1 — 2(s/t)qt = v^ 1 — 2q s for 
s G {0,t}, we get that Vs G [0,*], det(u _1 - 2(s/t)q t ) ^ <*=^> Vs G [0, i], det^" 1 - 2q s ) ^ 0, and thus 
25x,Q,6,Q;t n Sd(R) = 2\a:t H ^d(R)- Let v G T>x,a,b,a;t- Since I>a:,a,b,o;t is an open set, there is £ > such that 
u ± e7 d G f> x<atb<a . <t H Sd(R). Since 2? &!a;t is convex, u = (u + e7 d + w - e7 d )/2 G X> 6 ^ Q; i. □ 
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C Proofs of Section H 

C . 1 Proof of Theorem flOl 

Proof : From (|8]), we get: 

Li = ad{i ti y + 2x {i)i }9| i)i} + 2 ^ £{;,m}%,m} %,•>} + 2 X! ^{m,*}"^*,™}^} ( 46 ) 

l<m<d KmKii 

We want to show that LjL^ = LjLi for i 7^ j. Up to a permutation of the coordinates, Li and Lj are the 
same operators as L\ and L2. It is therefore sufficient to check that L±L 2 = L 2 L 2 . By a straightforward but 
tedious calculation, we get 
L\Li = 



a 2 



<9{i,i}<9{ 2 ,2} + 2a>x {2t2} d {1A} df 2 2} + 2a>Y / x {2,j} d {iA} d {2,2}d{2,j} + 

— ^ j¥2 



(o) (i) 

(2) 

f (^{1,2} + x U,k} d {hi} d {2, 3 } d {2,k}) + 2aaj {i,i} 9 {i,i} 9 {2,2} + 4z {M} a; { 2,2}<9f 1 ,1}<9{2,2} 

Z ' jjt2,kji2 * * ' " • ' 

(3) * , ' (1) (5) 

(4) 

A Y X {l,l} X {2, 3 } d {l.,l} d {2,3} d {2,2} + X{l,l}( 2d {l,n d {l,2} + Y X UM 9 { 1,1} 9 {2 J} 9 {2,k} ) + 

J"#2 ' v ' &1,k£2 

* * ' (?) » , ' 



(6) (8) 
2a Z{l,TO}<9{l,l}<9{i im }<9{2,2} + 4 X! S {l,mF{2,2}<9{l,l}<9{l,m}<9{2,2} + 

S v ' V v ' 

(2) (6) 
4 ( Y X {l-/m}X{ 2 .j}d{ 1 . 1} d{ 1 . m} d{ 2 . J }d{ 2 . 2} + Z{l,2}d{l,l}#{l,2}d{2,2}) 

> , ' (10) 

(9) 

+ X] :r {l,m}a :: {j,/c}5{i i i}9{ lim }9{2J}5{2,A;} + ^ 3f{l,m}^{l,2}^{l,m} 
mjtl,kjt2,j^2 m^l,m/2 

S v ' S v ' 

(11) (12) 

+2 ^2 2; {l,m}5{l i l}9{l i 2}9{2,m} + X {1,2}^{1,2} + 2:E {l,2}5{l i l}9{ li 2}9{2,2} 

> . ' (14) 

(13) 

^ ^{m,i}3{l ) m}3{l,/}^{2,2} + ^ x {2.2} x {-m,l}d{l.rn}d{i.l}df 22 y + 
m^l,/^l m^l,i^l 



(4) (8) 

2J a5{2,j}2;{m,i}^{l,m}9{l,/}^{2,2}3{2,j} + 2 a ; {2,m}^{l,2}5{i, m }9{ 2 ,2} + ^{2.2)^\\ ,2}^{2,2} + 
m^l,l^l,jjt2 ro^l,m^2 v 



(7) 



(11) (13) 



J Yj X {m,l} X {j,k}d{l,m}d{l,l}d{ 2 ^ k }d{ 2 .j} + ^2 :E {m,i}^{l,i}^{l,2}5{2,m} + X {2,m}#{ 12 } #{2,m} • 



jV2,fe#2 m^2,2^2 ■> 



(15) (16) 



(12) 
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Obviously, we have the same formula for L2L1 simply by exchanging the index 1 and 2. It is then sufficient 
to check that the above formula remains unchanged when we exchange the two index. Each term above is 
marked with a number. If this number is in the form (r), then the associated term is symmetric. Otherwise, 
there exist in the formula its corresponding symmetric term which is marked with the same number. □ 



C.2 Proof of Proposition 1111 

Proof : Let Xf ~ WISd{x,a,0, 12;t). We will check that for any polynomial function / of the matrix 
elements, we have E[f(Xf)] = E[f(X"'"' )]. Let us consider a polynomial function / of degree m: 

x e S d (R), f(x) = ^2 a 7 x 7 , 

7gN <i ( d + 1 )/ 2 ,|7|<m 

where [7] = J2i<i<j<d l7{*,i}l anc ^ ^ 7 = Yii<i<j<d x {i 'fy ■ Since the operator are affine, it is easy to check 
that Lf(x) and Lif(x) are also polynomial functions of degree m. We set: 

||/||p= |a 7 | and |L| = max \\Lx^\\ Pl 

7eN d(d + l)/2 ; | 7 |< m ' C 'I'l- 

so that ||L fe /||p < |L| fc |j/|jp for any A; € IN. Therefore, the series J2kL t k L k f(x)/k\ converges absolutely. By 
using I + 1 times Ito's formula, we get: 

£[f(xf )} = J2 y Lk fW + f { -^rr^nL l+1 f(x:)]d s . 

k=Q ' ^° 

Wishart processes have bounded moments since the drift and diffusion coefficients have a sublinear growth. 
Thus, C = max 7eN d( ti +i)/2 ; | 7 | <m sup se [ t ] E[|Xf 7 |] < 00 and we obtain that | J* Q E[L l+1 f(X* )]ds\ < 

C\\f\\ P {t\L\) l+1 /(l + 1)1 4 0. Thus, we have E[/(Xf)] = J2T=o t k L k f(x)/k! and similarly we get that 

r n,... x t \i Y n-l,...t 



f{X?" )\xl 



+K n 1, 



k 

fe„=0 



Now, we remark that C = maXyg^dfd+ii/a , |-y| <m su P s g [o,t] rnax(E[|X t ' x |], . . . , E[|X"'"' * |]) < 00 by using 

x l,m 

once again that Wishart processes have bounded moments. Since E[|L^ Tl /(X™ -1 ' " * )|] < C , ||/||p|L„|' c ™, 
we can switch the expectation with the series and get (fl~4")) . Then, since L k l n f(x) are polynomial function of 
degree m, we can iterate this argument and finally get (fT5j) , which gives the result. □ 



C.3 Proof of Theorem [H 

The proof is divided into two parts. First, we prove that the SDE (|24j) has a unique strong solution which 
is given by (|25|) and is well defined on 5j~(R). Second, we show that its infinitesimal generator is equal to 
the operator L\ defined in (p~6|) . 
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First step. Let us assume that (Xf )t>o is a solution to We use the matrix decomposition of 

(%i,jh<i,j<d given by ([23]) and set: 

r 

{Ut) {u +i} = E^MP^fM+i}. I €{!,..., r}, 

i=l 

r / r \ ^ r 

(c/ t ) { i,i} = w){i,i> - E E^^M^Hd.m} - (*?){!,!} - E(( c/ *){i^i}) 2 - 

/=i \j=i / ;=i 

We get by using Lemma [35] that: 

(^){M}+ELl((^){l,fe+l}) 2 ((^*){l,I+l})l<I<r 




((Ct){l,J+l})l<I<r ^ I I C T r k T T 

o / V 7 d _ r _i 

^) {1 ,l}+S =1 ((P t ){l,W}) 2 ((^){l,m})f<K^ ((^) {1 ,W})f<Kr^ \ 

c r((^t){l,;+l})l<i<r C r C^ C r /C, T = Xf . 

^r((^t){l,/+l})l<;<r KcJ / 

10 \ 

Since [ c r ] is invertible, Xf £ 5j"(IR) if, and only if: 

/ (^){i,i}+ELi((^){i.,+i}) 2 Wt) { l,l}h<l<r + 1 0\ 
Vz€R d , T (W){(,l})2<Kr+l 4 2 (47) 

\ / 

r 

= zi 2 (^t){i,i} +E^+! + (^*){i,*fi}«i) 2 > 0. > o. 

1=1 

In particular, we get that (£7o){i,i} = > since x £ <Sj~(IR). Now, by Ito calculus, we get from ([24]) 

that d(C/*) { i, i+1} = E[=i ELi(^ 1 )i,('v),k< +1 = and 

d(*7 t ) { i,i } = (a - r)dt + 2^{U t ) {l , l} dW? + 2 E[ =1 ELi (^ 1 )i, fe (X) { i, fe+ i}^ +1 

-ELi2((^) { 1 , + i } W +1 

= (a-r)dt + 2 A /(J7 t ) {M} dW t 1 . 

Thus, the solution (Xf ) t > is necessarily the one given by ([23]) (pathwise uniqueness holds for ((£7").m ;})i<;< r +i> 
and especially for the CIR diffusion (U^s^n since a > d — 1 > r). Reciprocally, it is easy to check by Ito 
calculus that (J2SJ) solves ([2"3]), 

Second step. Now, we want to show that L\ is the infinitesimal operator associated to the process 
(Xf )t>o- It is sufficient to compare the drift and the quadratic covariation of the process Xf with L\. Since 
the drift part of (Xf )t>o clearly corresponds to the first order of L%, we study directly the quadratic part. 
From fl2U, we have for i, j £ {2, . . . , d} 2 : 

d((xf) {M} ,(xf) {M} ) = 4((xf) {M} -n =1 [ELi(^ 1 )M(^f) { i, i+ i}] 2 + ELiELi(^ 1 ) fe ,/(xf) {1 , J+1} ] : 

= 4(Xf) {M} di, 

d((Xf ) {1 . i} , (Xf = ELll^l.-U^)]-!,!:* = (cC T )i-i d -idt = (X?){ itj ydt, 

d((Xf) {1 . 1} ,(Xf) {M} ) = 2EL 1 ELi(^)-i.fe( c r- 1 )M(^f){i.i+i}* = 2(Xf) {1 . l} dt, if *<r+l, 

d((Xf) {1 , 1} ,(Xf) {M} ) = 2EUELl(fcr)i-l-r ) fe(^ 1 ) fc ,/(Xf) {1>m} df 

= 2 E[=i( fc rC r 1 )i-i_ r ,z(Xf = 2(Xf) {M} dt if i > r+ 1, by Lemma[35] 

Thus, we deduce that L\ is the infinitesimal generator of (Xf )*>o- 
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Lemma 35 — Lety G Sj(R). We setr = Rk((j/ i> j) 2 <i,j<d), y\ = (j/i,i+i)i<i<r and y r { d = (yi, l+ i) r +i<i<d- 
We assume that there are an invertible matrix c r and a matrix k r defined on M.d—r—ixrffi) > such that 



(yi,jh<i,j<d 



C r 
k r 







Then, 


7 r,d 

we nave y 1 = 




V- 












o \ 







Proof 


: We set p = 


(t 


c r and have p 1 = 





c~ v 








&r Id—r—l / 


Vo 





. Since the matrix 



p-'yip-T = 



2/1,1 




-krC-hl) T ' 




I r 





r.d i — 1 r 

Vl ~ K rC r Hi 





o 



is positive semidefinite, we necessarily have y\' d — k r c r Y y\ = 0. 



D Proofs of Section [3] 

Lemma 36 — Let us consider an operator L on D that satisfies the required assumption and such that b 
and a have a sublinear growth, i.e. 3K > 0, ||6(a;)||2 + ||cr(a;)||2 < if (1 + ||x||2)- Let us consider Xf a potential 
vth-order scheme for L, with v > 1. Then, condition (|31|) holds. 



Proof : Wc have D C R c , and we set for p G IN*, = (eLi 2 ^ = (IMl2) 2p , ^ G D. Clearly, 

•/p G C^^D), and since is a potential ^th-order scheme for L, we have 

lim [E[/ P (Xf)] - f p {x)]/t = Lf p (x). (48) 

i->0+ 

Since we have dif p (x) = 2pxif p -\(x) and didjf p (x) = 4p(p — l)xiXjf p -%[x) + l;=j (2p/ p -i(x)), it is easy to 
check that \Lf p (x)\ < Y,Ui + IWWIWI^ 1 + E^i ^ 2 K\\ + \\xhf\\x\\ 2 r 2 < C p (l + \\x\\ 2 2 p ) = 

C p (l + f p (x)), for some constant C p > 0. From (|48|). we get that: 

3r)>0,Vte(0,r)),E[f P (X?)}<f p (x)+tC p (l + f p (x)). 

We consider now the scheme (X^,,i = 0,...,N) with N > T/n that starts from xq G D. We have 
E[/ p (X^ +i )] < E[/ p (X^)] + C p (T/7V)(l + E[/ p (X t p]). Let u = f p {x ) and = Ui(l + C p T/N) + C p T/N. 
We have E[f p {X$, )} < m and Ui = (l + C P T/Ny(f p {x ) + 1) - 1 < e c » T (f p {x Q ) + 1). Therefore we have 

sup sup nW^Wf] <e c - T (||a; ||^ + l) < oo, 

N>T/r,0<i<N 



which gives ([31 
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D.l Proof of Theorem [231 

Theorem 1231 defines the scheme as tt t h r (U")ir . We can prove (see later) that TJf is a potential ^th- 
order scheme for some operator, and it would be then easy to analyze the error if h r were in C££,i(<Sj~(R))- 
Unfortunately, h r is only smooth w.r.t. to the coordinates W{i.i}, . . . and is not a priori in Cw,i(«Sj~(R)). 

In fact, this is sufficient to show that 7r T /i r (£/")7r is a potential ^th-order scheme because these coordinates 
are the only one that are modified by the scheme. This requires however some further analysis which is made 
below. 

Let D C R^ be a domain. We introduce for any domain D C R*, £ € N* the set 

C£°oi| D (DxD) = {/eC°°(BxD,R), V 7 G N c , 3C 7 > 0, e 7 e IN*, V(x,x) € D, \d 7 f(x,x)\ < C 7 (l+|| (x, i)|| e ^)}, 

where ||.|| is a norm on R^ + ^ and <9 7 = 9J 1 , . . . ,<9^ c denotes the derivatives w.r.t. the coordinates of D. For 
/ G (D x D), we will say that (C 7 ,e 7 ) 7SN c is a good sequence for / if it is such that \d 7 f(x,x)\ < 

C 7 (l + \\(x : x)\\ e ^) for any (x,x) eDxD. 

Let us now consider an operator L that satisfies the required assumption on D. It is easy to check that 
all the iterated functions L k f are well defined and belong to C^ ol \ (D x D) as soon as / G C~ 1 L(D x D). 
Let us fix x, £ D and consider X* sampled according to a potential weak z/th-order scheme for L. Since 
x i ^ f(x,x) belongs to C~j(D), we know by (j3"0"]) that there are constants Cx,Ex,r]x such that 

<C^ +1 (l + \\x\\ E '). 

In practice, one would like instead to get some bounds where the dependence with respect to x is more 
tractable. This is why we introduce the following definition. 



Vt G (0,7fc 



^ k\ 

k=0 



t k L k 



f(x,x) 



Definition 37 — Let L be an operator that satisfies the required assumption on D. We will say that a 
potential weak vth-order scheme for L satisfies the immersion property if for any D C R e , £ G N* and any 
function f G C^ )1 | D (D x D) with a good sequence (C 7 ,e 7 ) 7SN c, there exist positive constants C, E, and ij 
depending only on (C 7 ,e 7 ) 7gN <; such that 



Vte (0,77), 



E[f(Jtf,x)]-£~1*L k f(x,x) 



fc=0 



<ct^(i + \\(x,x)\n. 



In practice, most of the usual schemes satisfy this property. In fact, to prove that a scheme is a potential 
z^th-order scheme, it is common to use a Taylor expansion that gives generally at the same time the immersion 
property. We illustrate this for the exact scheme below. 



Proposition 38 — Let D C R c ,, b : D -> R c and a : D -> A4 C (R) such that \\b(x)\\ + \\a(x)\\ < C(l + ||a;||) for 

some C > 0, and assume that Lf(x) = J2i=i bi(x)dif(x) + 1 5Z< j=\{ aaT { x ))i,j®i®jf{ x ) satisfies the required 
assumption on D. Then, for any v G IN, f/ie exaci scheme is a potential weak vth-order scheme for L and it 
satisfies the immersion property. 

Proof : Let / G C^,jL(D x D). We know from the sublinear growth condition that we have bounds on the 
moments oiXf: Vg G N*,3C, > 0,Vi G [0, 1], E[||(Xf,z)|| 9 ] < C q (l + \\(x,x)p). By iterating Ito's Formula, 
we get then easily for t G [0,1], 

E[f(Xf , x)} = £ ^L k f(x, x) + f ^Eli^+VK, Z)\d8. 
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Since L u+1 f g C^, 1 | D (D x D), there are constants C > and g £ W* depending only on a good sequence 
of / such that \L u+1 f(x,x)\ < C(l + ||(a;,f)||«). Thus, we deduce that \E[f(X?,x)] - ELo T\L k f(x, x)\ < 
^C(l + C 9 {l + \\(x,xW)). □ 

Besides, the immersion property is easily preserved by scheme composition. 

Proposition 39 — Let L\ and L2 be two operators that satisfy the required assumptions on D, and assume 
that p x (t)(dz) and p x (t)(dz) are respectively potential weak vth-order discretization schemes on D for these 
operators that satisfy the immersion property. Let Ai,A2 > and X x ° t '£ it ~ p 2 (\2t) o p x (X 1 t)(dz). Let 

f £ C^ )1 | D (D x D). Then, there are constants C,E,v that only depend on a good sequence of f such that 



Vie (0,7?), 



E ^t h+h ^L l if(x,x) 



<ct v+L (i+\\(x,x)\n. 



Therefore, 



If L1L2 = L2L1, p 2 (t) o p l x (t)(dz) is a potential weak vth-order discretization scheme for L\ + L2 
satisfying the immersion property. 

If v > 2, p 2 (t/2) op 1 (i) op 2 .(t/2) and \ (p 2 (t) o p x (t) +p 1 (t) op^(t)) are potential weak second order 
schemes for L\ + L2 satisfying the immersion property. 



This proposition is a straightforward extension of Proposition 1.15, Corollary 1.16 and Theorem 1.17 of 
Alfonsi PQ, and we do not repeat the proof here. Thanks to this result, we can prove the immersion property 
of the schemes that are obtained by splitting. The Ninomiya-Victoir scheme [20] which is obtained by a 
composition of exact schemes naturally satisfies this property. By looking at the proof of Theorem 1.18 
in PP, it still satisfies this property if we replace the Gaussian samples by moment matching variables. Also, 
we can check that the second and third order schemes for the CIR process presented in Alfonsi [1] satisfy 
the immersion property. 

Corollary 40 — The Ninomiya- Victoir scheme and the second and third order scheme for the Cox-Ingersoll- 
Ross process given in satisfy the immersion property. 



Corollary 41 — Let L\ (resp. L2) be an operator that satisfies the required assumptions on Di (resp. 
D2J. Let Xl ,Xl and Xt' X2 be potential weak vth-order schemes for L\ and L2 sampled independently. Then, 
(X]' Xl , Xt' X2 ) is a potential weak vth-order schemes on Di x D2 that satisfies the immersion property. 

Proof : From the immersion property, it is easy to check that (Xl' Xl , X2) (resp. {x\, Xf ,X2 )) is a potential 
^th order scheme for L\ (resp. L2) on Di x D2 that satisfies the immersion property. The composition of 
these schemes is simply (X t ' X1 , X^' X2 ). Since L\ and L2 operate on different domains, we have L1L2 = L2L1, 
which gives the result. □ 



Proof of Theorem [ 

Let x € S+(R) and / e C~ ol (<S+(R)) and r = Rk((x hJ ;) 2 <i,j<d)- Since the operator L\ satisfies the required 
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assumption, we know that the exact scheme is a potential i/th-order scheme (Proposition |3"8")) , and there are 
constants C, E, n > depending on a good sequence of / such that 

Vt G (0, i,), |E[/(Xf)] - £ yLlf(x)\ < C(l + (49) 

fc=0 

On the other hand, we know from Theorem [TBI equation (f28|) and Corollary [15] that we have 

X x t = n T h r (U?)n, 

where (J/")/! n solves the SDEs (f26|) starting from the initial condition ity for 1 < I < r + 1, and (£^™){i,j} = 
uuj\ for the other coordinates. We have also Xf = 7r T h r (Ut)ir by construction, and it is natural to focus 
on the function u n- f(ir T h r (u)ir). 
Let us consider the set 

{x e S d (R), s.t. {Xi tj )2<i,j<d e SjljjRj.xi,! > 0}. 

It is isomorphic to (R+ x R d_1 ) x 5jl x (R) by the map x i-> • ■ ■ , {xi,j)2<i,j<d)- We have to notice 

now that the function h r defined by ([25]) is such that /i r G C^ >1 | R+x[Rti _ 1 (R + x x <Sj"_ 1 (R)). It is indeed 

a polynomial function with respect to Mix, . . . , Then, it is easy to check that it i— > f(w T h r (u)ir) G 

C^il,, xKd _!(R + x R d_1 x 5jl x (R)) since / G C™ ol {s£(R)). Moreover, by the chain rule, we can get a good 
sequence for this function that only depend on a good sequence of / since ir and h r are fixed. 

By assumption, (£/"){!,!} is a potential j/th-order scheme for the operator (a — r)9{i : i} + 2u{ lil }9^ 1 ^ 

and satisfy the immersion property. The schemes (2 < i < r + 1) can be seen as a Ninomiya- 

Victoir scheme with moment matching variables. They are therefore potential z/th-order scheme for the 
operator ~ ([I], Theorem 1.18) and satisfy the immersion property from Corollary 1401 Therefore, from 

Corollary HU ((C r "){i,i}, • • ■ , iPt){i,d}) 1S a potential uth order scheme for (a — r)d{i t iy + 2u{i ) i}9| 1 u + 
\ Yli—i i} satisfying the immersion property. Thus, there are constants that we still denote by C, E, r\ > 
depending on a good sequence of / such that: 

vt g (0,7?), |E[/(7r T M(y t »] - q/(7r T M[) t >)]| < + IMI £ )- (50) 

Now, one has to notice that ||u|| < C (1+|| a; ||) for some constant C" > since we have U{i,i}+2fe=i( u {i,k+i}) 2 — 
and x and x have the same Frobenius norm. We get then the result by gathering (j49|) and (fBTJl) . □ 



Remark 42 — W^e can check by looking at the proof above that the scheme obtained for L\ also satisfies the 
immersion property. Since we do not need it for the construction of the further schemes for WISd(x, a, b, a) 
or AFFd(x,a, B,a), we will no longer mention it. 

D.2 Proof of the third condition of Theorem 1211 for Wishart processes 

The scope of this part is to show the following result. 



Proposition 43 — Let (Xf) t > ~ WIS d (x,a,b,a), f G C™ ol (S d (R)), x G 5+(R) and T > 0. The 
u(t,x) = E{f(Xt)] is C°° on [0,T] x 6>j~(R) ; solves dtu(t,x) = Lu(t,x) and its derivatives satisfy 



VI G N,Vn G N^^aC ; ,„,e/,„ > 0,Vx G 5+(R), Vt G [0,T], 



n 



l<i<j<c 



< C L n(l+\\x\\ e '-). 



(51) 
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For technical reasons on the boundary of <Sj~(R), it is convenient here to work with / G C^ >1 (5 £ ;(IR)) rather 
than / G C™ ol {Sj(R)). The proof of Proposition 1351 is based on a remarkable identity on the derivatives of 
the Laplace transform which is stated below. 

Lemma 44 — Let (X*) t >g ~ WISd{x,a,b,a) and v — Vr + ivj such that vr G 2\ ;t and Vj G Sd(R)- 

— Law 

We denote by 4>(t,a,x,v) the Laplace transform of Xf given by (jlOp , the other parameters a, b being fixed. 
Then, the derivative w.r.t Xs^n satisfies the following equality 

®{k,i}<i>(t,a,x,v) = <j>(t,a + 2,x,v)p\ k ^ (v), 

where p\ k ^ is a polynomial function of the matrix elements of degree d defined by : 



Pt ( v ) = Tr «adj(7 d - 2q t v)m t {e k / + t k ^ie d ,k )mf =: ^ a]* k ' l *v y , where v 1 = vjl'jf ■ 

7GN 2 ,\f\<d X Jt 

Moreover, its coefficients are bounded uniformly in time: 

(K.w>|)<*«. 



k.l 



l,k\ T 



3K t > 0, Vs G [0,f 

Proof : We get from (fTU|) , 
d{k,i}<j>{t,a,x,v) 



max 

d(d+l) 

7GN 2 , |t| <d 



Tr 



wadj(/ d - 2q t v)m t {e d + t k ^ie d ' )mj 



exp (Tr [i>(J d - 2q t v) 1 m t xmf] ) 



det(I d - 2q t v) 



det(/ d - 2g t u) 2 



= <f)(t, a + 2,x, v)Ti v&d]{I d - 2q t v)m t {e k / + lk^ie l ^)mi 

Since s i— > ||m s |j and s n- ||g s || arc continuous functions on [0,t], we obtain the bounds on the polynomial 
coefficients. □ 



l,k\ T 



Proof of Proposition \%Bf Let / G C^ ol (Sd{R))- First, let us observe that (|5T|) is obvious when I = \n\ = 
Since we have VZ G M,L l f G C^i(S d (R)), and d l t u{t,x) = £(L l f(Xf)), it is sufficient to prove (J5lj) only for 
the derivatives w.r.t. x. 

We first focus on the case \n\ = 1 and want to show that d[k,i}u(t,x) satisfies ff5T|) - The sketch of 
this proof is to write / as the inverse Fourier transform of its Fourier transform and then use Lemma 1441 
Unfortunately, / has not a priori the required integrability to do that, and we have to introduce an auxiliary 
function f p . 

Definition of the new function /„. Since 2\ a; T given by © is an open set and G 2?& Tj there 
is p > such that pld G l\ a ;T- Let fx : R — > R be the function such that [i(x) = if x < — 1 or 
x > 0, fi(x) — exp( a ,^ 1 +1 ^ ) if — 1 < x < 0. We have \i G C°°(R). We consider then the cutoff function 

C : R -> R G C°°(R) defined as Vx G R, C(z) = ^f^y^y ■ U is nondecreasing, such that < ((a;) < 1, 
C(x) — if x < — I and ((x) = 1 if x > 0. Besides, we have C G C£^i(R) since all its derivatives have a 
compact support. Now, we define a # G C^ >1 (5 ( j(R)) as 

i? : S d (R) -> R 
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It is important to notice that < i? < 1, = 1 if x G 5d(R) and t?(x) = if there is i € {1, . . . , d} 

such that < —1 or i < j G {l,...,d} such that x 2 ^-^ > 1 + x^^jx^^j . Let 7 e ug^- 1 )/ 2 . Since 

/ € ^(^(R)), there are constants if, J5 > and if', £?' > such that, Vx € S d (R) 

|<9 7 (t9/)(:z:)| < if(l + ||x|| B )]J (l{|x {M} |>-l} 

) n 

l<i<j<d 

< ^ , (i+n(!c {<lt - } )i< i < <J ii^)nt 1 ( i {i-{M>i>-i}) n (^,. 1} <i+« { M>«u.,>}) • 

l<i<j<<2 

Here, the upper bound only involves the diagonal coefficients. We define 

x G 5 d (R), / p (x) := t?(x)/(x)exp(-Tr(px)), 

and obtain from the last inequality that f p belongs to the Schwartz space of rapidly decreasing functions 
since p > 0. Thus, its Fourier transform also belongs to the Schwartz space and we have 

fp( x ) = ~ — , d ( d+i) / d ( d+i) exp(-Tr(wx))J r (/ p )(i;)dw, where T{f p )(v) = I d(d+1) exp(Tr(ivx))f p (x)dx, 

(27r) 2 JR 2 JR 2 

and in particular f p ,F{f p ) G L x («S d (R)) n L°°(5 d (R)). 

A new representation of ii(t,x). We have /(x) = exp(pTr(x))/ p (x) for x G Sj~ (R), and therefore 

tx(t,x) = E[exp(Tr(pXf))/ p (Xf)] 
1 



d ( d+i) 
[zir) 2 



E 



exp(Tr[(- l « + P I d )Xf])F{f p ){v)dv 

R 2 



. E[exp(Tr[(-z« + p/ d )Xf])]^(/ p )(z;)dz;. 

(27TJ 2 JR 2 

The last equality holds since J d(d+p |E[exp(Tr[(— it; + pI d )Xf ])] | x | J r (/ p )(t;)|dt; < cp(t, a, x, pId)\\J-(f p )\\i < 

\R 2 

00. Here, we have used that pld G £>b, a ;T to get <fi(t, a, x, pld) < 00. 

Derivation with respect to X{k.i} 5 k,l G {1, . . . , c£}. From Lemma 1441 we have by Lebesgue's theorem 

d{ k d}u{t,x) = - — 1 j cj)(t,a + 2,x,-iv + pI d )p\ kA (pld ~ iv)F(f p )(v)dv (52) 

(2tt)— — Jr 2-^ 

since \df kl} (/>(t, a, x, -iv + pI d )T(f p )(v)\ < \4>(t,a + 2,x, pl d )\\p\ k,l} (pl d - iv)T(f p )(v)\ and p\ k ' l} [pl d - 
iv)J-(f p )(v) is a rapidly decreasing function. 

Let 1 < k', I' < d. An integration by part gives J R (pId — iv){k',i'} exp (Tr[x(iv — pld)]) $(x)f(x)dx{k',i'} — 
+ l fc/=i ,) f R ex-p(Tr[x(iv - pl d )])d {k , d , } {i9(x)f(x))dx {k , d , } , and thus 

(p/ d -it;) { ^,// } ^(exp[-pTr(x)]t?(x)/(x))(t;) = + l^ =i 0J r (exp[-pTr(x)]9 {v , n [t?(x)/(x)])(t;). 

We set ^(7) = rii</ ! '<r<d(^ ii + lfc'=(') 7{fc '' i,} for 7 G N d(d+1)/2 and get by iterating the argument that: 
[] (ph - iv?$$F{f p ){v) = ^(7)J-(exp[-pTr(x)] 9 7 (t9 x /)(x)) (v). (53) 

l<fe'</'<d 

Since p\ k ' 1 ^ (pl d - iv) = dja+n a 7,{fc,Z} IIi<fc'<z'<<i(W<i ~~ ^Jl'V? ' we S et from and ([53]): 

7SN 2 , 1 T I < <i — — — 1 > J 

d {k<l} u(t,x) = Yl ^' {fc "V(7)E(a 7 (/ x 0)(Y t *)) = Y, «7' {fc ''V(7)E(a 7 /(^)) , (54) 
h\<d h\<d 
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where (Yf)t>o ~ WISd(x,a + 2,6, a). Here, we have used that <9 7 (i? x f)(y) = d 1 f(y) for y e <?j~(IR). 

Law 

From Lemma 23] (a?'^' ) d(d+i) is bounded for t € [0, T], and we get ([51]) when Inl = 1 since 

7GN 2 ,|7|<d 

d-yf G C^,j(»Sd(R)). Thanks to (j54"l) . a derivative of order |n|, can be seen as a (bounded) linear combination 
of derivatives of order \n\ — 1, and we easily get (I51[) by an induction on \n\. 

It remains to check that we have indeed dtu(t,x) = Lu(t,x). Let t, h > 0. By the Markov property, we 
have w(i + h, x) ~ E[u(t, X%)]. From (f5~Tj) and Ito's formula, we get [u(i + ft., a;) — u(t, x)]/h — > a;). □ 



Lemma 45 — Lei a, x G Sj~(R), B G £(<S+(R)) that satisfies ©, and a;(i) oe i/ie solution of the following 
ODE 

x(t)=x + / (a + B(a;(s)))(te. (55) 
j o 

Then, we have x(t) E S^(R) for t > 0. 

Proof : The ODE ([55]) is afhne and has unique solution on Sj~(R) which is given by 

t > 0, x(t) = exp(tB)(x) + [ exp(sB) (a)ds, (56) 

Jo 

where V£ G R+ Vx € S d {R), exp(tB)(x) = J2T=a » Bfe (- T ) = 5 ° • • • such that B °( a; ) = x - 

/c times 

We assume first that a,x G <Sj~'*(R) and consider r = inf{i > 0, x(£) ^ 5j~(R)}, with the convention 
inf = +oo. We have r > 0. Let us assume by a way of contradiction that r < oo. Then, x(t) cannot be 
invertible and there is y G 5j~(R) such that y ^ and Tr(yx(r)) = 0. From (f56|) and (j4}, we get 

Tr{x'{T)y) = Tr ([B(x(r)) + a]y) > 0, 

since a is positive definite. Therefore, there is e G (0, r) such that Tr(yx(r — e)) < 0. Let us recall now that 
z G <Sj~(R) <*=^> Vy G 5j"(R), Tr(yz) > 0. Thus, x(t - e) £ Sj~(R), which contradicts the definition of r. 

In the general case a, a; € <5>j~(R), we observe that the solution (|56|) is continuous w.r.t. a; and a, and 
thus yt > 0, x(t) G Sj~(R) since Sj~(R) is a closed set. □ 
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